陳銘，張士誠，胥云，馬新仿，鄒雨時(shí)

（1.中國石油大學(xué)（北京），北京 102249；2.Texas A & M University，College Station 77840，USA；3.中國石油勘探開發(fā)研究院，北京 100083；4.中國石油油氣藏改造重點(diǎn)實(shí)驗(yàn)室，河北廊坊 065007）

0 引言

水平井分段多簇體積改造是非常規(guī)油氣開發(fā)的關(guān)鍵技術(shù)[1-3]。為提高儲量動用程度，北美油公司不斷探索縮小井距與縫間距的體積改造技術(shù)[4]?，F(xiàn)場光纖溫度監(jiān)測與聲監(jiān)測表明[5]，小間距改造情況下，多簇裂縫存在不均衡進(jìn)液等問題，同時(shí)各簇進(jìn)液比例與射孔參數(shù)、各簇地應(yīng)力分布等因素緊密相關(guān)。為提高多簇裂縫均衡擴(kuò)展程度、優(yōu)化多簇壓裂方案設(shè)計(jì)，急需建立一種耦合“井筒-射孔-裂縫”的多裂縫擴(kuò)展高效模擬方法[6]。

裂縫擴(kuò)展模型包括二維、擬三維、平面三維和全三維模型[6]。二維模型以PKN（Perkins-Kern-Nordgren）和KGD（Khristianovich-Geertsma-Daneshy）模型為代表，適用于恒定縫高的單縫擴(kuò)展。擬三維模型包括橢圓模型與基于PKN單元的模型。橢圓模型假設(shè)裂縫由以射孔為中心的上下半橢圓構(gòu)成，裂縫形態(tài)僅需要由橢圓長軸、短軸和離心率確定[7]?；赑KN單元的擬三維模型引入縫高解析解計(jì)算裂縫高度，流動仍為沿縫長的一維流動，計(jì)算量較小，但在射孔段為高應(yīng)力、薄互層等情況下縫高誤差較大[8]。Kresse等[9]提出的非常規(guī)裂縫模型即為擬三維模型在多裂縫擴(kuò)展中的應(yīng)用。趙金洲等[10]提出了基于擬三維模型的多裂縫擴(kuò)展模型，并建立了射孔優(yōu)化方法。擬三維模型在多簇應(yīng)力干擾計(jì)算方面仍是一種二維模型方法，難以對應(yīng)力干擾下縫高擴(kuò)展進(jìn)行準(zhǔn)確模擬[11]。為準(zhǔn)確分析裂縫形態(tài)，Advani[12]、Barree等[13]提出平面三維模型。平面三維模型采用三維固體方程計(jì)算巖石變形，通過裂縫邊界確定裂縫長度和高度。Peirce等[14]引入尖端解析解，提出了平面三維模型的隱式水平集算法。Dontsov等[15]提出尖端統(tǒng)一解析解，并應(yīng)用于隱式水平集算法[16]。Chen等[17]采用隱式水平集算法對比了拉鏈?zhǔn)綁毫雅c同步壓裂的裂縫擴(kuò)展形態(tài)。為描述水力裂縫空間扭轉(zhuǎn)問題，Carter等[18]提出了全三維模型，全三維模型可計(jì)算裂縫空間扭轉(zhuǎn)，計(jì)算量巨大。Xu等[19]實(shí)現(xiàn)了平面偏轉(zhuǎn)的分段多簇壓裂模擬，并研發(fā)了FrackOptima軟件。全三維模型目前工業(yè)應(yīng)用尚未見報(bào)道，主要困難是理論上空間扭轉(zhuǎn)裂縫的判斷準(zhǔn)則還不完善[20]，同時(shí)計(jì)算量巨大，不利于工程應(yīng)用。此外，近幾年研究者也提出了天然裂縫發(fā)育地層的復(fù)雜裂縫擴(kuò)展模型[21-24]。本文重點(diǎn)研究多簇壓裂中“井筒-射孔-水力裂縫”耦合流動的多裂縫擴(kuò)展問題，并分析多簇壓裂的射孔設(shè)計(jì)對策，暫不考慮天然裂縫的影響[10,14,16,25]。

顯然，目前在分段多簇壓裂設(shè)計(jì)中最為實(shí)用與可靠的壓裂模型為平面三維模型。目前廣受認(rèn)可的算法為 Peirce等[14]、Dontsov等[16]隱式水平集算法。該算法采用隱式方法求解多裂縫擴(kuò)展的流固耦合方程，裂縫邊界通過隱式水平集方法判斷。隱式解法雖然可以保證計(jì)算穩(wěn)定，但非線性方程組的求解工作量較大，尤其是邊界元系數(shù)矩陣為稠密矩陣，更加增大了計(jì)算量，同時(shí)在處理支撐劑運(yùn)移、裂縫閉合等問題時(shí)方程約束增多，進(jìn)一步增大了求解難度和計(jì)算量，不利于工程應(yīng)用。此外，該算法目前并未考慮射孔摩阻作用，相關(guān)工程應(yīng)用報(bào)道也較少[25]。

為此，本文提出一種平面三維模型求解新算法，該算法采用顯式積分方法求解裂縫擴(kuò)展的流固耦合方程，通過最短路徑算法與尖端解析解求解裂縫邊界。通過與解析解[26]、澳大利亞國家研究院有機(jī)玻璃實(shí)驗(yàn)[27]、隱式水平集算法[16,28]進(jìn)行對比，驗(yàn)證算法的準(zhǔn)確性與效率。采用浙江油田頁巖氣水平井實(shí)際參數(shù)進(jìn)行算例分析，重點(diǎn)分析簇間應(yīng)力分布和射孔數(shù)分布對各簇進(jìn)液控制規(guī)律及工程對策。

1 數(shù)學(xué)模型

1.1 基本假設(shè)

本文研究多簇壓裂工藝的裂縫擴(kuò)展規(guī)律，幾何模型如圖1所示，每段簇?cái)?shù)為nf，裂縫簇序號用k標(biāo)記。多簇裂縫擴(kuò)展主要包含 4個(gè)物理過程：流體在井筒和裂縫內(nèi)流動；巖石在流體壓力作用下變形；流體向地層濾失；裂縫尖端破裂。本文采用 Peirce等[14]平面三維模型的假設(shè)，即濾失符合 Carter濾失模型；縫內(nèi)流體邊緣與裂縫邊緣重合；裂縫沿垂直于最小水平主應(yīng)力方向的平面擴(kuò)展。Bunger等[29-31]研究表明，多裂縫擴(kuò)展中偏轉(zhuǎn)并不明顯，同時(shí)國內(nèi)非常規(guī)油氣儲集層水平應(yīng)力差較大，因此平面模型適用。

圖1 幾何模型

1.2 控制方程

1.2.1 固體方程

用固體方程描述流體壓力、地應(yīng)力作用下裂縫寬度的分布。根據(jù)無限大地層彈性力學(xué)點(diǎn)源解，固體變形的邊界積分方程為[14]：

1.2.2 流動方程

①井筒流動方程。井筒內(nèi)流動摩阻對流量分配影響較小[32]，同時(shí)各簇間距較小，因此各簇之間井筒流動摩阻可以忽略。各簇裂縫對應(yīng)的井底壓力滿足：

射孔摩阻計(jì)算公式為[33]：

根據(jù)流量守恒，各簇流量之和為總注入排量：

②縫內(nèi)流動方程。對于每簇裂縫，縫內(nèi)流體遵循泊肅葉流動，即：

考慮流體不可壓縮，因此流體連續(xù)性方程為：

將（5）式代入（6）式，得到流體流動控制方程：

1.2.3 裂縫邊界條件

根據(jù)斷裂力學(xué)準(zhǔn)則，裂縫尖端達(dá)到擴(kuò)展條件時(shí)：

由于裂縫尖端即為流體邊緣，因此裂縫邊界滿足零流量條件，即：

1.2.4 裂縫尖端解析解

斷裂力學(xué)解僅限于裂縫尖端很小的范圍，因此需要非常細(xì)的網(wǎng)格才能準(zhǔn)確捕捉（8）式所示的尖端條件。為實(shí)現(xiàn)粗網(wǎng)格條件下對尖端條件的準(zhǔn)確捕捉，本文采用縫尖解析解求解裂縫位置[16]。圖2為尖端斷裂力學(xué)解與解析解的適用范圍示意圖。

圖2 尖端解析解與斷裂力學(xué)解適用范圍示意圖

引入無因次量，令：

裂縫尖端擴(kuò)展速度、擴(kuò)展長度、斷裂韌性、彈性模量、濾失系數(shù)等均影響裂縫擴(kuò)展的臨界寬度。裂縫尖端擴(kuò)展時(shí)滿足方程：

（11）式為裂縫擴(kuò)展速度、擴(kuò)展步長與臨界寬度的非線性方程。由于δ的變化范圍較小，為簡化方程求解，采用 Dontsov提出的近似解法，近似解誤差在0.3%以內(nèi)[15]。具體解法如下。

①取δ=0時(shí)，得到（11）式的零階近似：

②根據(jù)（12）式計(jì)算δ：

③將（13）式的計(jì)算結(jié)果代入（11）式，得到修正的零階近似解：

（14）式為尖端速度、臨界寬度、擴(kuò)展步長的控制方程，該方程適用范圍遠(yuǎn)大于斷裂力學(xué)解范圍，因此可有效增大空間步長，提高計(jì)算效率。

2 求解算法

2.1 網(wǎng)格系統(tǒng)

本文采用固定網(wǎng)格計(jì)算裂縫擴(kuò)展。將裂縫平面劃分為足夠多的矩形單元，單元類型包含 4種：注入點(diǎn)單元、已開啟單元、縫尖單元和未開啟單元（見圖3）。每個(gè)時(shí)間步需判斷縫尖單元是否達(dá)到擴(kuò)展條件，從而更新網(wǎng)格的單元類型。單元中心點(diǎn)為寬度和壓力求解點(diǎn)，單元邊界為流量求解位置。

圖3 網(wǎng)格系統(tǒng)與單元類型

2.2 離散方程

2.2.1 固體方程的離散

固體方程離散采用常位移不連續(xù)法[34]。單元中心點(diǎn)寬度近似為單元寬度，因此邊界積分方程的離散形式為

（15）式的矩陣形式為：

2.2.2 流動方程的離散

①縫內(nèi)流動方程的離散。對控制方程（7）在空間內(nèi)進(jìn)行離散，得到流動方程的一階微分方程形式：

（17）式中，A(w)p與S的分量形式分別為：

②井筒流動方程的離散。井筒流動方程的離散形式為：

上式為nf+1維非線性方程組，F(xiàn)的分量形式為：

2.3 流固耦合方程

微分方程（17）中θ=0時(shí)為顯式格式，θ=1時(shí)為完全隱式格式。微分方程（17）為剛性方程，顯式方法需滿足CFL（Courant-Friedrichs-Lewy）條件才能保證計(jì)算穩(wěn)定，因此研究者主要采用隱式方法求解微分方程（17）[6-7]。隱式解法滿足無條件穩(wěn)定性，不需要CFL條件，但隱式方法在流固耦合方程計(jì)算中仍然存在問題，主要表現(xiàn)在：①求解方程（17）通常采用牛頓-拉夫遜或不動點(diǎn)迭代方法[6-7]，每次迭代需要解線性方程組。由于固體方程中的影響系數(shù)矩陣C為稠密矩陣，因此通常采用直接法解線性方程組。直接法的工作量取決于單元數(shù)量，單元數(shù)量增大后，時(shí)間復(fù)雜度巨大。盡管預(yù)處理方法可以將工作量降低，但對于單元數(shù)量較大情況，時(shí)間復(fù)雜度仍然巨大[35]。②考慮支撐劑運(yùn)移或裂縫閉合等非線性問題時(shí)，流固耦合方程的約束條件增多，增大了隱式方法求解難度和計(jì)算工作量[35]。③為保證計(jì)算精度，隱式方法不能采用太大的時(shí)間步長[7,28]。

考慮到隱式方法仍存在計(jì)算困難和時(shí)間成本較高等問題，本文采用顯式方法求解方程（17）。為提高顯式方法計(jì)算效率，采用穩(wěn)定型 RKL（Runge-Kutta-Legendre）[36]顯式積分方法求解方程（17）。該方法采用多步時(shí)間積分方法，利用遞歸式Legendre多項(xiàng)式的穩(wěn)定性特征，擴(kuò)大了顯式算法的CFL時(shí)間步長，因此可顯著降低計(jì)算量，提高顯式算法計(jì)算效率。

方程（17）可記為：

利用Legendre遞歸多項(xiàng)式絕對值小于1的特征，將單步積分采用S步遞歸式積分求解。時(shí)間項(xiàng)2階精度的S步 RKL方法的計(jì)算格式參照文獻(xiàn)[36]確定。S步RKL方法的時(shí)間步長滿足：

對于裂縫擴(kuò)展問題，在保證穩(wěn)定性條件下，確定時(shí)間步長時(shí)也需考慮計(jì)算精度。因此取S步RKL方法的時(shí)間步長為：

（24）式中ε為松弛系數(shù)，0＜ε≤1，經(jīng)過計(jì)算分析，其取0.8可滿足足夠精度。對于非均質(zhì)地應(yīng)力問題，若裂縫擴(kuò)展進(jìn)入低應(yīng)力層后擴(kuò)展速度加快，需要適度降低松弛因子以避免時(shí)間步過大。

2.4 尖端擴(kuò)展

Peirce等[14]、Dontsov等[16]利用與裂縫尖端相鄰的激活單元的解析解計(jì)算裂縫尖端位置和尖端寬度。本文顯式算法以尖端未開啟單元臨界寬度和擴(kuò)展速度計(jì)算單元開啟狀態(tài)。每一時(shí)間步判斷縫尖單元的相鄰未開啟單元是否達(dá)到擴(kuò)展條件，從而更新單元類型。采用最短時(shí)間路徑算法確定單元開啟時(shí)間與是否達(dá)到開啟條件。需要注意的是，本文仍采用網(wǎng)格激活方式進(jìn)行裂縫邊界捕捉[13]。

以圖4為例，對于待開啟單元a，相鄰單元為b、c、d、e。根據(jù)最短路徑算法，a單元開啟時(shí)間為：

未開啟單元的當(dāng)前開啟時(shí)間為正無窮大。假設(shè)b、c、d單元為未開啟單元，e為已開啟單元，則a單元開啟時(shí)間為：

當(dāng)vea滿足尖端擴(kuò)展條件（11）式時(shí)，a單元成為開啟單元。

圖4 裂縫尖端擴(kuò)展示意圖

2.5 算法

平面三維多裂縫擴(kuò)展模型算法如下。

①設(shè)定注入時(shí)間tf，模型輸入?yún)?shù)包括：巖石力學(xué)參數(shù)、注入程序、液體參數(shù)、地應(yīng)力分布等；

②令α=0，t=t0，計(jì)算初始寬度w0和壓力p0；

③令α=α+1，t=t+Δt，求解（22）式得到當(dāng)前開啟單元的寬度wα和壓力pα。若t＞tf，結(jié)束計(jì)算；

④采用牛頓-拉夫遜方法求解（20）式得到各簇流量分配；

⑤計(jì)算當(dāng)前時(shí)刻的待開啟單元數(shù)，求解（11）式得到所有待定單元的臨界寬度，判斷是否發(fā)生開啟，若發(fā)生開啟，則更新為尖端單元，否則仍為未開啟單元；

⑥根據(jù)步驟⑤更新單元類型，確定新的裂縫尖端單元和待定單元，返回步驟③。

算法主要包括各簇分流量計(jì)算、固體變形與流動耦合方程計(jì)算和裂縫擴(kuò)展邊界計(jì)算 3個(gè)模塊。各簇分流量通過牛頓-拉夫遜方法計(jì)算，具有二階收斂速度。固體變形與流動耦合方程的穩(wěn)定型RKL解法具有二階時(shí)間精度，并滿足計(jì)算穩(wěn)定性[36]。裂縫擴(kuò)展邊界采用尖端解析解計(jì)算，尖端解析解適用范圍可達(dá)縫長的10%，而線彈性力學(xué)解析解范圍僅為縫長的 0.1%～1.0%，因此可在粗網(wǎng)格條件下實(shí)現(xiàn)較高計(jì)算精度，從而減小計(jì)算量[37]。3個(gè)模塊均滿足收斂性和穩(wěn)定性，因此算法理論上可行。

3 算法驗(yàn)證與效率對比

3.1 準(zhǔn)確性驗(yàn)證

3.1.1 與penny裂縫解析解對比

為驗(yàn)證算法準(zhǔn)確性，首先與解析解對比?？紤]無層間應(yīng)力差情況，該情況下單裂縫擴(kuò)展符合penny模型。礦場條件下，penny裂縫擴(kuò)展主要為黏性主導(dǎo)能量耗散方式，因此采用黏性主導(dǎo)penny裂縫進(jìn)行單縫驗(yàn)證。

單縫驗(yàn)證的基本參數(shù)為：排量 5 m3/min，流體黏度5 mPa·s，彈性模量30 GPa，泊松比0.2，斷裂韌性0.2 MPa·m1/2，濾失系數(shù)0，注入時(shí)間10 min。徑向裂縫擴(kuò)展的特征時(shí)間為[30]：

當(dāng)注入時(shí)間遠(yuǎn)小于特征時(shí)間時(shí)，裂縫擴(kuò)展能量消耗以縫內(nèi)流動摩阻耗散為主；當(dāng)注入時(shí)間遠(yuǎn)大于特征時(shí)間時(shí)，裂縫擴(kuò)展能量消耗以尖端破裂能量耗散為主，兩者中間為過渡過程。計(jì)算特征時(shí)間為3.30×109min，則注入時(shí)間（10 min）遠(yuǎn)小于特征時(shí)間，為黏性主導(dǎo)裂縫。黏性主導(dǎo)徑向裂縫半徑和裂縫入口寬度為[30]：

該算例采用單元尺寸為2.5 m×2.5 m，計(jì)算結(jié)果如圖5所示。結(jié)果顯示，本文算法計(jì)算結(jié)果與解析解吻合，表明算法可準(zhǔn)確計(jì)算penny裂縫擴(kuò)展動態(tài)。

3.1.2 與隱式水平集算法對比

為進(jìn)一步驗(yàn)證算法準(zhǔn)確性，考慮分層加載應(yīng)力、流體發(fā)生濾失的情況，采用 Dontsov等[16]隱式水平集算法進(jìn)行驗(yàn)證。設(shè)置3層最小主應(yīng)力，如圖6a所示。其他參數(shù)[16]為：彈性模量9.5 GPa，泊松比0.2，流體黏度 0.1 Pa·s，注入排量 0.01 m3/s，斷裂韌性 1 MPa·m1/2，濾失系數(shù) 2.065×10-6m/s1/2。

圖5 本文算法計(jì)算結(jié)果與penny裂縫解析解對比

圖6 分層加載應(yīng)力、流體發(fā)生濾失情況下的裂縫擴(kuò)展剖面與Dontsov等[16]隱式水平集算法裂縫擴(kuò)展輪廓對比

圖6b為注入3 600 s時(shí)本文算法得到的裂縫寬度剖面與文獻(xiàn)[16]裂縫輪廓對比。結(jié)果顯示，本文結(jié)果與Dontsov等[16]隱式水平集算法結(jié)果吻合，表明本文算法可以準(zhǔn)確求解存在濾失的裂縫擴(kuò)展形態(tài)。

3.2 計(jì)算效率對比

采用算例為澳大利亞國家科學(xué)院有機(jī)玻璃實(shí)驗(yàn)[27]。有機(jī)玻璃滿足均質(zhì)線彈性，并可以拍照監(jiān)測裂縫擴(kuò)展動態(tài)。具體實(shí)驗(yàn)參數(shù)為：流體黏度 30 Pa·s，有機(jī)玻璃彈性模量3.3 GPa，泊松比0.4。分3層加載水平應(yīng)力，如圖7所示。實(shí)驗(yàn)采用變排量注入[27]。圖8為實(shí)驗(yàn)665 s時(shí)的裂縫形態(tài)與本文算法計(jì)算結(jié)果的對比?？梢钥闯?，本文算法計(jì)算結(jié)果接近實(shí)驗(yàn)獲得的裂縫形態(tài)，進(jìn)一步驗(yàn)證了本文算法準(zhǔn)確性，也說明可以采用實(shí)驗(yàn)參數(shù)進(jìn)行計(jì)算效率對比。

圖7 Wu等[27]有機(jī)玻璃壓裂實(shí)驗(yàn)應(yīng)力加載與裂縫形態(tài)

圖8 本文算法與澳大利亞研究院實(shí)驗(yàn)結(jié)果對比

Zia等[28]采用隱式水平集算法及文獻(xiàn)[27]中實(shí)驗(yàn)參數(shù)進(jìn)行了計(jì)算，并評價(jià)了算法計(jì)算效率。本文采用與文獻(xiàn)[28]相同大小網(wǎng)格（0.43 cm×0.43 cm），處理器均為Intel(R) Core(TM) i7-5600 CPU@2.60GHz。本文采用MATLAB編程，文獻(xiàn)[28]采用Python編程，兩種語言在科學(xué)計(jì)算方面執(zhí)行效率相當(dāng)，因此通過對比兩者CPU占用時(shí)間確定本文算法與隱式水平集算法的計(jì)算效率。本文算法采用25、50與75級積分算法占用CPU時(shí)間均在300 s以內(nèi)，其中75級積分算法占用CPU時(shí)間僅為73 s，而文獻(xiàn)[28]隱式水平集算法占用CPU時(shí)間為619 s?？梢姳疚乃惴ㄓ?jì)算效率大幅提高。

4 實(shí)際井算例分析

以浙江油田昭通頁巖氣示范區(qū)YS112H4水平井組YS112H4-1井為例進(jìn)行實(shí)例分析。YS112H4-1井目的層為下志留統(tǒng)龍馬溪組，采用小間距分段多簇壓裂工藝，設(shè)計(jì)平均簇間距10 m，施工排量14 m3/min，施工液體為FAB-2滑溜水，密度為1 016 kg/m3，黏度為2 mPa·s。儲集層厚度為34 m，水平應(yīng)力差為17.0～19.5 MPa。巖石彈性模量為 35.7 GPa，泊松比為 0.27，斷裂韌性為1.0 MPa·m1/2。儲集層基質(zhì)滲透率 1.0×10-7μm2，因此可以忽略流體向基質(zhì)濾失。目標(biāo)層段上下有一定應(yīng)力遮擋，水平最小主應(yīng)力剖面如圖9所示。圖9中z=0 m為射孔位置所在深度?？籽壑睆?2 mm，射孔修正系數(shù)0.7。

圖9 地應(yīng)力剖面與分簇射孔示意圖

由于已有較多文獻(xiàn)[38-40]分析分簇?cái)?shù)、射孔直徑等參數(shù)與各簇進(jìn)液量的關(guān)系，本文重點(diǎn)分析實(shí)例井簇間應(yīng)力分布、射孔數(shù)及射孔數(shù)分布等對各簇進(jìn)液量和裂縫擴(kuò)展的影響。

4.1 簇間應(yīng)力分布對各簇進(jìn)液的控制作用

4.1.1 簇間應(yīng)力均質(zhì)分布

首先分析地應(yīng)力平面均質(zhì)的情況，即各射孔簇最小水平主應(yīng)力相同，均為60 MPa。該情況下射孔摩阻與縫間應(yīng)力干擾控制各簇進(jìn)液分布。

由圖10可知，射孔簇應(yīng)力均質(zhì)分布情況下，每簇4、8、16孔的射孔方式各簇均能開啟并進(jìn)液。每簇4、8孔時(shí)，各簇進(jìn)液量接近；每簇16孔時(shí)，外側(cè)簇裂縫進(jìn)液為中間簇裂縫進(jìn)液的1.24倍。因此，單簇孔數(shù)越多，各簇進(jìn)液量差異越大。

由圖11可知，每簇4、8孔時(shí)，各簇裂縫長度接近；每簇16孔時(shí)，外側(cè)簇裂縫半長為330 m，而中間簇裂縫半長為265 m，前者為后者的1.24倍，與進(jìn)液量比例一致。由于產(chǎn)層上部應(yīng)力小于產(chǎn)層下部應(yīng)力，裂縫高度傾向于向產(chǎn)層上部擴(kuò)展，而產(chǎn)層以下基本沒有裂縫展布。中間簇近井區(qū)域裂縫高度大于外側(cè)簇，說明近井區(qū)域的應(yīng)力干擾最為顯著，近井區(qū)域中間裂縫長度擴(kuò)展受阻，因而容易發(fā)生縱向擴(kuò)展。

圖11 簇間應(yīng)力均質(zhì)分布時(shí)的裂縫擴(kuò)展形態(tài)

4.1.2 簇間應(yīng)力非均質(zhì)分布

設(shè)定射孔簇1—5所在層位的最小水平主應(yīng)力分別為62，60，60，60，60 MPa，即射孔簇1最小水平主應(yīng)力比其他射孔簇高2 MPa。

由圖12可知，每簇4或8孔時(shí)，高應(yīng)力簇進(jìn)液量低于其他簇，但各簇均能進(jìn)液；每簇16孔時(shí)，高應(yīng)力簇不能開啟進(jìn)液。說明每簇16孔不能平衡簇間2 MPa應(yīng)力差，每簇8孔可以實(shí)現(xiàn)各簇有效進(jìn)液，每簇4孔各簇進(jìn)液幾乎沒有差別。

結(jié)合簇間應(yīng)力均質(zhì)分布的進(jìn)液情況可知，盡管射孔簇 1裂縫為應(yīng)力干擾作用下的進(jìn)液主導(dǎo)縫，但若其最小水平主應(yīng)力高于其他簇 2 MPa，則可能無法開啟進(jìn)液。計(jì)算結(jié)果表明，簇間應(yīng)力非均質(zhì)對各簇進(jìn)液的控制作用大于應(yīng)力干擾。因此，簇間應(yīng)力分布是分簇限流設(shè)計(jì)的重要參數(shù)。

圖12 簇間應(yīng)力非均質(zhì)分布時(shí)的各簇進(jìn)液量分布

圖13 簇間應(yīng)力非均質(zhì)分布時(shí)的裂縫擴(kuò)展形態(tài)

由圖13可知，每簇4、8孔時(shí)盡管各簇進(jìn)液量差別不顯著，但由于產(chǎn)層上部地應(yīng)力小于產(chǎn)層下部地應(yīng)力，同時(shí)裂縫 1長度擴(kuò)展受到其他裂縫應(yīng)力干擾“擠壓”作用，因此裂縫 1傾向于向產(chǎn)層上部擴(kuò)展，進(jìn)而出現(xiàn)縫高過量擴(kuò)展。裂縫擴(kuò)展形態(tài)受層間應(yīng)力剖面和應(yīng)力干擾共同控制。水力裂縫會選擇阻力最小路徑擴(kuò)展，當(dāng)簇間應(yīng)力干擾作用大于層間應(yīng)力差作用時(shí)，裂縫會選擇縱向擴(kuò)展路徑。受層間應(yīng)力分布和應(yīng)力干擾影響，高應(yīng)力射孔簇的裂縫形態(tài)發(fā)生改變，進(jìn)而影響改造效果。這表明各簇產(chǎn)量與進(jìn)液量的關(guān)系不一定一致。這種現(xiàn)象也難以通過井下光纖溫度或聲監(jiān)測識別，因此在設(shè)計(jì)階段應(yīng)對簇間應(yīng)力分布進(jìn)行細(xì)致解釋分析，并盡量將應(yīng)力接近的區(qū)域劃分為一段。

4.2 射孔數(shù)分布對各簇進(jìn)液的控制作用

4.2.1 簇間應(yīng)力均質(zhì)分布

Cramer等[41]通過現(xiàn)場試驗(yàn)井下拍照發(fā)現(xiàn)，由于射孔質(zhì)量、射孔相位差別，壓裂過程并非每個(gè)射孔均可開啟，因此有必要分析不均勻射孔數(shù)對各簇進(jìn)液的影響。設(shè)計(jì) 3種射孔方案進(jìn)行對比：①各簇射孔數(shù)依次為8，9，8，8，8；②各簇射孔數(shù)依次為8，10，8，8，8；③各簇射孔數(shù)依次為8，11，8，8，8。設(shè)定各射孔簇最小水平主應(yīng)力相同，均為60 MPa。

由圖14可知，射孔簇2增加1～3孔之后，該簇進(jìn)液最多。增加孔數(shù)越多，該簇進(jìn)液量占比越大。增加1，2，3個(gè)孔時(shí)，射孔簇2進(jìn)液量占總注入量的比例分別為21.5%，23.3%和25.0%。增加3個(gè)射孔后，射孔簇2與其他各簇進(jìn)液量差異較大，因此應(yīng)避免增加3個(gè)以上射孔數(shù)。由圖15可知，某簇射孔數(shù)增加過多，會導(dǎo)致該簇裂縫縫長過大，簇間改造反而更不均衡?？梢?，增大某一簇的射孔數(shù)會提高該簇的進(jìn)液量，從而平衡地應(yīng)力、近井摩阻或應(yīng)力干擾差異的影響，但射孔數(shù)增加不應(yīng)過多，增加1～2孔的作用已經(jīng)非常顯著。本文模擬結(jié)果與Cramer等[41]現(xiàn)場試驗(yàn)得到的認(rèn)識一致。

圖14 簇間應(yīng)力均質(zhì)分布時(shí)不同射孔方案的各簇進(jìn)液分布

4.2.2 簇間應(yīng)力非均質(zhì)分布

設(shè)定射孔簇1—5所在層位的最小水平主應(yīng)力分別為62，60，60，60，60 MPa，即射孔簇1最小水平主應(yīng)力比其他射孔簇高2 MPa。為增加射孔簇1的進(jìn)液比例，設(shè)計(jì)3種射孔方案：①各簇射孔數(shù)依次為9，8，8，8，8；②各簇射孔數(shù)依次為10，8，8，8，8；③各簇射孔數(shù)依次為11，8，8，8，8。

圖15 簇間應(yīng)力均質(zhì)分布時(shí)不同射孔方案的裂縫擴(kuò)展形態(tài)

圖16 簇間應(yīng)力非均質(zhì)分布時(shí)不同射孔方案的各簇進(jìn)液分布

圖17 簇間應(yīng)力非均質(zhì)分布時(shí)不同射孔方案的裂縫擴(kuò)展形態(tài)

由圖16可知，盡管射孔簇1為高應(yīng)力簇，但增加該簇1～3個(gè)射孔后，該簇進(jìn)液量與其他4簇差異逐漸減小，增加該簇 2個(gè)射孔即可實(shí)現(xiàn)均勻進(jìn)液。但從圖17可以看出分析，高應(yīng)力簇的改造仍然不充分，裂縫1縫高出現(xiàn)過量擴(kuò)展。由于裂縫1在產(chǎn)層內(nèi)擴(kuò)展受到其他簇較大應(yīng)力干擾作用，水力裂縫會沿最小阻力路徑擴(kuò)展，因此裂縫1傾向于向產(chǎn)層上部擴(kuò)展，導(dǎo)致裂縫1縫高擴(kuò)展過大，產(chǎn)層改造面積不足。同時(shí)，增加 3孔與增加 2孔的裂縫形態(tài)較為接近。研究表明，增加高應(yīng)力簇2個(gè)射孔數(shù)可以提高該簇進(jìn)液，達(dá)到均勻進(jìn)液；均勻進(jìn)液并不意味著裂縫形態(tài)均勻，裂縫形態(tài)受應(yīng)力干擾和層間應(yīng)力剖面兩方面控制。體積改造多簇壓裂設(shè)計(jì)應(yīng)加強(qiáng)對層間應(yīng)力剖面、層內(nèi)平面應(yīng)力分布的解釋，并應(yīng)盡量將應(yīng)力接近的地層作為一段進(jìn)行分簇壓裂。

5 結(jié)論

提出了水平井分段壓裂平面三維多裂縫擴(kuò)展模型求解新算法，算法準(zhǔn)確可靠，與目前廣受認(rèn)可的隱式水平集算法相比，新算法計(jì)算速度大幅提高。以浙江油田昭通頁巖氣示范區(qū)下志留統(tǒng)龍馬溪組頁巖氣水平井實(shí)際參數(shù)進(jìn)行模擬分析，研究發(fā)現(xiàn)：各簇射孔數(shù)相同時(shí)，簇間應(yīng)力非均質(zhì)對各簇進(jìn)液的控制作用大于應(yīng)力干擾，減小單簇射孔數(shù)可以平衡簇間應(yīng)力差異，實(shí)現(xiàn)均衡進(jìn)液；調(diào)整各簇射孔數(shù)量可以實(shí)現(xiàn)均衡進(jìn)液，各簇射孔數(shù)差別應(yīng)控制為1～2孔，簇間射孔數(shù)差別太大，會引起液體過多進(jìn)入射孔數(shù)最多的簇，反而加重各簇進(jìn)液不均勻，不利于各簇均衡改造；增加高應(yīng)力簇的射孔數(shù)有利于均勻進(jìn)液，但進(jìn)液均勻并不等于裂縫形態(tài)一致，裂縫形態(tài)受應(yīng)力干擾和層間應(yīng)力剖面共同控制；水力裂縫沿最小阻力路徑擴(kuò)展，對于高應(yīng)力射孔簇，若層間應(yīng)力差較小，高應(yīng)力簇裂縫更容易沿縱向擴(kuò)展，不利于產(chǎn)層內(nèi)儲集層改造，體積改造應(yīng)選擇應(yīng)力接近的區(qū)域作為一段。

致謝：感謝國家留學(xué)基金委的資助，感謝中國石油勘探開發(fā)研究院王臻博士在算法設(shè)計(jì)方面給予的幫助。

符號注釋：

A——裂縫覆蓋區(qū)域面積，m2；A——流動方程系數(shù)矩陣；A(t)——t時(shí)刻的裂縫覆蓋區(qū)域面積，m2；Cg——邊界積分方程格林函數(shù)，具體形式可參見文獻(xiàn)[34]，Pa/m3；C——影響系數(shù)矩陣，Pa/m；Clv——濾失系數(shù)，m/s1/2；dk——第k簇裂縫的射孔直徑，m；E——巖石彈性模量，Pa；G——剪切模量，Pa；h，r——網(wǎng)格單元在z方向的序號；i，m——網(wǎng)格單元在x方向的序號；i0，j0，h0——注入點(diǎn)單元在x，y，z方向的序號；j，n——網(wǎng)格單元在y方向的序號；k——裂縫簇序號；K——射孔磨蝕修正系數(shù)；KIc——Ⅰ（張）型裂縫斷裂韌性，Pa·m1/2；M(w)——（17）式等號右端項(xiàng)對應(yīng)的系數(shù)矩陣；nf——每段簇?cái)?shù)；nk——第k簇裂縫的射孔數(shù)；p——縫內(nèi)流體壓力，Pa；p——縫內(nèi)流體壓力矩陣，Pa；p0——初始時(shí)刻的縫內(nèi)流體壓力，Pa；p0——初始時(shí)刻的縫內(nèi)流體壓力矩陣，Pa；pin,k——第k簇裂縫的裂縫入口壓力，Pa；pp,k——第k簇裂縫的射孔摩阻，Pa；pα——tα?xí)r刻的縫內(nèi)流體壓力，Pa；pw——井底壓力，Pa；q——單位長度體積流量矢量，m2/s；Qk——第k簇裂縫的進(jìn)液流量，m3/s；Qt——注入排量，m3/s；R——裂縫半徑，m；s——距縫尖的距離，m；S——源匯項(xiàng)系數(shù)矩陣；t——時(shí)間，s；t0(x,y,z) ——坐標(biāo)(x,y,z)處發(fā)生濾失的時(shí)刻，s；t0a，t0b，t0c，t0d，t0e——a、b、c、d和e單元開啟時(shí)間，s；tc——徑向裂縫擴(kuò)展的特征時(shí)間，s；tf——注入時(shí)間，s；tα——第α?xí)r間步，s；Δt——時(shí)間步長，s；ΔtE——顯式歐拉差分格式時(shí)間步長，s；v——裂縫尖端擴(kuò)展速度，m/s；vba，vca，vda，vea——b、c、d和e單元到a單元的擴(kuò)展速度，m/s；w——裂縫寬度，m；w——裂縫寬度矩陣，m；w0——初始時(shí)刻的裂縫寬度；w0——初始時(shí)刻的裂縫寬度矩陣；win——裂縫入口寬度，m；wα——tα?xí)r刻的裂縫寬度，m；wα——tα?xí)r刻的裂縫寬度矩陣，m；——裂縫平均寬度，m；x，y，z——場點(diǎn)坐標(biāo)，m；x'，y'，z'——源點(diǎn)坐標(biāo)，m；Δx，Δy——空間步長，m；α——時(shí)間步編號；δ——常量，無因次；δk(x,y,z)——第k簇裂縫的狄拉克函數(shù)，m-2；θ——常系數(shù)，0≤θ≤1；μ——流體黏度，Pa·s；ρ——液體密度，kg/m3；σh——最小水平主應(yīng)力，Pa；σh——最小水平主應(yīng)力矩陣，Pa；Δτ——RKL方法的時(shí)間步長，s；υ——巖石泊松比。