基于能量守恒的半解析半數(shù)值多裂縫擴(kuò)展模擬方法

鐘安海, 魯明晶*, 張潦源, 錢欽, 賀文君, 程呈, 朱海燕

(1.勝利油田分公司石油工程技術(shù)研究院, 東營(yíng) 257000; 2.成都理工大學(xué)能源學(xué)院, 成都 610059)

水平井分段多簇壓裂是開(kāi)發(fā)油氣資源的有效增產(chǎn)措施之一,研究壓裂過(guò)程中水力裂縫的擴(kuò)展延伸對(duì)油氣的高效開(kāi)采是至關(guān)重要的[1],這直接關(guān)系到油氣井的產(chǎn)量。

為此,中外學(xué)者對(duì)多簇裂縫的競(jìng)爭(zhēng)擴(kuò)展進(jìn)行了大量研究,總結(jié)出多種方法。曹楷楠等[2]結(jié)合近場(chǎng)動(dòng)力學(xué)和有限差分法模擬了多裂縫擴(kuò)展。程萬(wàn)等[3]采用邊界元法建立了水平井多簇裂縫同步擴(kuò)展模型,研究了多裂縫的競(jìng)爭(zhēng)擴(kuò)展機(jī)制。郭建春等[4]采用二維位移不連續(xù)法建立了誘導(dǎo)應(yīng)力場(chǎng)分布數(shù)學(xué)模型對(duì)多裂縫擴(kuò)展進(jìn)行了研究。朱海燕等[5-6]基于離散裂縫方法模擬研究了頁(yè)巖油藏與氣藏的多裂縫擴(kuò)展情況。許江文等[7]和焦戰(zhàn)等[8]通過(guò)擴(kuò)展有限元方法模擬了裂縫擴(kuò)展。Dontsov等[9]提出了隱式水平集算法能夠捕捉裂縫擴(kuò)展的多尺度行為并定位動(dòng)態(tài)移動(dòng)的裂縫前緣。Cheng等[10-11]采用非線性節(jié)理單元建立流固耦合多裂縫閉合模型,對(duì)水平井分段壓裂的裂縫擴(kuò)展形態(tài)進(jìn)行了模擬。Li等[12]擴(kuò)展了基于滲透率的水平集方法(permeability-based hydraulic fracture, level set method,PHF-LSM),模擬了非均質(zhì)巖層中的多簇裂縫擴(kuò)展。Huang等[13]基于有限元-離散裂縫網(wǎng)絡(luò)模型研究了工程因素對(duì)多裂縫垂向擴(kuò)展的影響。Hunsweck等[14]采用有限元法,Olson等[15]采用擴(kuò)展有限元法模擬了水平井中多裂縫同步擴(kuò)展問(wèn)題。Wu等[16]基于位移不連續(xù)法建立二維多裂縫擴(kuò)展模型,模擬了水平井多簇裂縫的競(jìng)爭(zhēng)擴(kuò)展問(wèn)題。上述常見(jiàn)的裂縫擴(kuò)展數(shù)值模擬方法,都在不同程度上存在模擬精度不高或者運(yùn)算效率較低的缺點(diǎn),且多局限于二維尺度,難以獲取水力裂縫整體形態(tài),在模擬工程尺度下的多裂縫競(jìng)爭(zhēng)擴(kuò)展問(wèn)題上具有一定的局限性。

為此,提出一種將無(wú)量綱近似解與能量守恒相結(jié)合的裂縫擴(kuò)展模擬方法。該方法可充分考慮水力壓裂過(guò)程中的濾失影響及多裂縫間的競(jìng)爭(zhēng)擴(kuò)展效應(yīng)。在此基礎(chǔ)上,編制半解析半數(shù)值裂縫擴(kuò)展程序C5Frac,并將計(jì)算結(jié)果與ILSA II[17-18]的模擬結(jié)果進(jìn)行對(duì)比驗(yàn)證。該方法可以實(shí)現(xiàn)多裂縫擴(kuò)展的快速模擬,且具有極高的計(jì)算精度,能夠?yàn)閴毫逊桨冈O(shè)計(jì)提供實(shí)時(shí)性指導(dǎo)。

1 半解析半數(shù)值裂縫擴(kuò)展模型

1.1 控制方程

如圖1所示,水力裂縫垂直于水平井筒向前擴(kuò)展。模型為單個(gè)壓裂段,段長(zhǎng)為Z,共包含N條裂縫,兩兩裂縫間的距離分別為hk(k=1,2,…,N-1),則可得出關(guān)系式為

圖1 單個(gè)壓裂段內(nèi)多裂縫擴(kuò)展示意圖Fig.1 Schematic of multiple fractures propagation in a fracture stage

(1)

裂縫面的徑向擴(kuò)展由井筒處注入的不可壓縮牛頓流體所驅(qū)動(dòng),并在具有滲透性的線彈性巖石中準(zhǔn)靜態(tài)延伸(即遠(yuǎn)低于巖石聲速),其特征為

(2)

式(2)中:E為楊氏模量;ν為泊松比;E′為平面應(yīng)變假設(shè)下的巖石特征參數(shù)。

韌性K′與巖石的斷裂韌性KIC的關(guān)系可表示為

(3)

為簡(jiǎn)化上述問(wèn)題,還需引入以下假設(shè):①裂縫擴(kuò)展遵循線彈性斷裂力學(xué);②所有裂縫面都呈徑向延伸且相互平行;③重力在彈性方程和流體流動(dòng)方程中均被忽略;④流體前緣與裂縫尖端重合,即不考慮流體滯后問(wèn)題;⑤遠(yuǎn)場(chǎng)地應(yīng)力均勻且恒定;⑥不考慮裂縫的扭曲和轉(zhuǎn)向。

提出上述假設(shè)之后,對(duì)于每條裂縫,需要求解的未知量包括:裂縫寬度w(r,t)、縫內(nèi)壓力pf(r,t)、縫內(nèi)體積流量q(r,t)、裂縫半徑R(t)、其他裂縫施加的作用力σ(r,t)以及縫口流量Q(t)。其中,r為裂縫內(nèi)某一點(diǎn)到裂縫中心的距離,t為時(shí)間?？刂品匠倘缦?。

(1)采用卡特濾失定律描述縫內(nèi)流體向地層的濾失,則縫內(nèi)流體的連續(xù)性方程為

(4)

式(4)中:CL為濾失系數(shù);k為巖石滲透性;cr為儲(chǔ)層壓縮系數(shù);φ為巖石孔隙度;μ為流體動(dòng)力黏度;σo為最小水平主應(yīng)力;po為儲(chǔ)層孔隙壓力;t0為水力裂縫某處發(fā)生濾失的時(shí)刻。

(2)縫內(nèi)流動(dòng)方程由泊肅葉流動(dòng)方程描述為

(5)

(3)通過(guò)非局部積分關(guān)系將裂縫寬度w(r,t)與牽引力T相互耦合得到裂縫的彈性方程為

(6)

式(6)中:x和s均為裂縫內(nèi)某一點(diǎn)到裂縫中心的距離;T為裂縫受到的牽引力,由縫內(nèi)部壓力、其他裂縫施加的相互作用力和遠(yuǎn)場(chǎng)地應(yīng)力組成,其表達(dá)式為

T(ρ,t)=pf(ρR,t)-σ(ρR,t)-σo

(7)

(4)根據(jù)線彈性斷裂力學(xué),裂縫的Ⅰ型應(yīng)力強(qiáng)度因子為[19]

(8)

當(dāng)應(yīng)力強(qiáng)度因子達(dá)到Ⅰ型斷裂韌性KIC時(shí),裂縫即向前擴(kuò)展。

(5)根據(jù)每條水力裂縫內(nèi)部壓力的分布,將其他裂縫施加在裂縫i上的壓應(yīng)力進(jìn)行疊加,計(jì)算裂縫i受到的總互作用力,計(jì)算公式為

(9)

式(9)中:σi為裂縫i受到的其他裂縫疊加的總互作用力;σj,i為裂縫j對(duì)裂縫i的互作用力。

(6)假設(shè)井筒內(nèi)流體流動(dòng)時(shí)的壓力損失為零,且流量之和等于總排量Qo,即井筒滿足體積平衡。因此得

(10)

式(10)中:Rw為井筒入口半徑;pf(N)為第N條裂縫的縫內(nèi)壓力。

因此,每條裂縫共包含:4個(gè)場(chǎng)方程、1個(gè)移動(dòng)邊界方程(裂縫擴(kuò)展條件)和1個(gè)流體守恒方程。

整個(gè)系統(tǒng)的初始條件為

R=0,w=0,q=0,pf=0

(11)

裂縫尖端的邊界條件為[20-21]

w(R,t)=0,q(R,t)=0

(12)

入口處邊界條件為

2πrq(r,t)=Q

(13)

在初始和邊界條件下可得到全局質(zhì)量平衡方程為

(14)

此外,還能得到雷諾潤(rùn)滑方程為

(15)

通過(guò)上述方程聯(lián)立求解,確定w(r,t)、pf(r,t)、q(r,t)、R(t)、σ(r,t)和Q(t)。

1.2 近似解

1.2.1 壓力分布近似解

函數(shù)形式同水力裂縫入口端和前緣附近的預(yù)期壓力漸進(jìn)性保持一致。壓力分布函數(shù)的假設(shè),通過(guò)消除基于流體流動(dòng)控制方程的有限差分離散化來(lái)計(jì)算每個(gè)時(shí)間步上的分布需要,可大大降低計(jì)算量。通過(guò)如下方程來(lái)表達(dá)流體壓力。

(16)

式(16)中:μc為所提出的“復(fù)合黏度”,該物理量反映濾失與韌性相關(guān)的額外能量損耗,是一個(gè)作用類似黏度的復(fù)合耗散參數(shù);Π(ρ,t)為無(wú)量綱壓力分布函數(shù),可表示為

(17)

式(17)中:A、B和w為系數(shù),在恒定泵注排量和無(wú)限均勻彈性巖石的假設(shè)下,分別取0.358 1、0.092 69和2.479;ψ(t)為空間均勻壓力,與韌性相關(guān)[22]。

1.2.2 韌性近似解

裂縫擴(kuò)展條件為KI=KIC,KI由牽引力T(ρ,t)計(jì)算得

(18)

(19)

因此,式(19)提供在任何時(shí)間都滿足隱式擴(kuò)展的壓力分布,使得復(fù)合黏度能夠在每一個(gè)時(shí)間步中進(jìn)行求解。并且,雖然增加了由裂縫擴(kuò)展給出的新控制方程,但通過(guò)顯式求解ψ(t),強(qiáng)制μc(t)與韌性關(guān)聯(lián),因此參與迭代的變量的數(shù)量沒(méi)有發(fā)生變化。

1.2.3 相互作用力近似解

裂縫延伸時(shí)的不均勻性和瞬態(tài)壓力的全彈性解是耗費(fèi)計(jì)算時(shí)間的主要原因。為實(shí)現(xiàn)高效計(jì)算,Cheng等[23]提出一種近似方法,其中非均勻壓力被均勻壓力所替代,在每個(gè)時(shí)間步中為每個(gè)水力裂縫選擇該均勻壓力,從而產(chǎn)生與非均勻內(nèi)壓擴(kuò)展出的實(shí)際水力裂縫體積相同的模擬裂縫,即

(20)

式(20)中:Pj為第j條均壓裂縫引起的內(nèi)部均勻靜壓力;wj為裂縫j的寬度。

相互作用力模型通過(guò)對(duì)所有相鄰裂縫的互作用應(yīng)力求和來(lái)實(shí)現(xiàn)。因此,施加在水力裂縫上的互作用力近似表示為

(21)

1.2.4 彈性近似解

局部裂縫寬度w(ρ,t)出現(xiàn)在體積平衡中,其中還包括入口邊界條件中使用的入口寬度。牽引力T(ρ,t)可由式(22)給出。

(22)

w(ρ,t)由半徑r決定,將半徑表示為無(wú)量綱半徑γ(t)的乘積以簡(jiǎn)化運(yùn)算。半徑通過(guò)求解近似方程組和特征半徑得

(23)

通過(guò)引入非局部彈性關(guān)系計(jì)算w(ρ,t)的所有必要變量,求解入口w(0,t)的開(kāi)度。將壓力代入泊肅葉方程中,得到另一個(gè)由泊肅葉定律導(dǎo)出的裂縫寬度。反過(guò)來(lái),由入口邊界條件得出約束限制為

(24)

1.2.5 全局體積平衡

通過(guò)對(duì)受初始和裂縫尖端邊界條件影響的局部體積平衡進(jìn)行積分,得到整體體積平衡方程為

(25)

式(25)中:C′L=2CL。

(26)

1.2.6 縫口條件

如式(10)所示,縫口條件由縫口壓力等值以及井筒流量守恒給出。滿足這些條件需要求出縫口壓力的近似值,這將從計(jì)算得出的流體壓力分布中獲取。對(duì)流體壓力分布進(jìn)行近似估算,由于函數(shù)形式在縫口處具有奇點(diǎn),計(jì)算縫口壓力需要限制井筒半徑,因此,解對(duì)井筒半徑的敏感性往往會(huì)由于縫口附近的壓力梯度較大而產(chǎn)生很大的誤差。因此將這些入口壓力視為未知數(shù),定義它們以符合整體能量守恒。

1.3 全局能量守恒

由于井筒壓力變化,很難準(zhǔn)確估計(jì)井筒壓力。因此,通過(guò)相反的近似來(lái)計(jì)算井筒的壓力,以滿足總能量守恒,將與韌性相關(guān)的能量添加到能量守恒式中,這使得新模型對(duì)所有的體系都有效,包括所謂的“韌性主導(dǎo)”體系。更新后的功率平衡為

(27)

式(27)中:pf(Rw,t)Q(t)為裂縫的能量輸入速率;Dc為與巖石壓裂有關(guān)的耗散率,可表示為

(28)

Df為與黏性流體流動(dòng)相關(guān)的耗散率,可表示為

(29)

DL為與濾失相關(guān)的流體損失率,可表示為:

(30)

U為通過(guò)變形巖石應(yīng)變能來(lái)增加應(yīng)變能的一部分,是可恢復(fù)的彈性能,可表示為

(31)

(32)

Wo為原地應(yīng)力對(duì)裂縫的作用。為考慮流體損失,原位應(yīng)力的功被修正為

(33)

式(33)中:S為用于量化非均勻的原位應(yīng)力;α為濾失相關(guān)系數(shù)。

Pperf為通過(guò)縫口的功率損失,可通過(guò)經(jīng)典壓降方程[式(34)]計(jì)算。

(34)

通過(guò)解壓力、寬度和半徑對(duì)流量的隱性依賴表達(dá)式代替未知的μc和γ得

(35)

(36)

(37)

(38)

(39)

(40)

1.4 模型求解流程

求得近似解的步驟如下。

步驟1定義輸入?yún)?shù)σc、CL、E′、μ′、KI、Qo、h。

步驟3使用牛頓法求解式(24)和式(25)。

步驟4求解能量守恒方程。

步驟6將時(shí)間步長(zhǎng)向前推進(jìn)Δt,Q(t)和α(t)可作為Q(t+Δt)和α(t+Δt)的預(yù)估值。

步驟7重復(fù)步驟3～步驟6,直至達(dá)到所需的總泵送時(shí)間。

2 模型驗(yàn)證及計(jì)算效率對(duì)比

為驗(yàn)證所建立的半解析半數(shù)值方法,將該方法模擬得到的計(jì)算結(jié)果與參考解結(jié)果進(jìn)行對(duì)比。根據(jù)參考解所具有的適用性,采用兩類不同的參考解分別對(duì)單裂縫擴(kuò)展和多裂縫擴(kuò)展進(jìn)行驗(yàn)證,并對(duì)比多裂縫擴(kuò)展模擬時(shí)的計(jì)算效率。

2.1 單裂縫擴(kuò)展

單裂縫擴(kuò)展的參考解是由Dontsov[24]開(kāi)發(fā)的一種數(shù)值近似解計(jì)算得出的。該方法可捕捉水力裂縫在低濾失與高濾失間以及高黏度與強(qiáng)韌性狀態(tài)間的過(guò)渡行為,還可捕捉以韌性或黏度為主的水力裂縫和以濾失和存儲(chǔ)為主的水力裂縫的擴(kuò)展行為,是一種可考慮濾失的裂縫擴(kuò)展高精度模擬方法。

圖2為C5Frac模擬結(jié)果與參考解的對(duì)比結(jié)果,可以看出,兩者的誤差在1%以內(nèi),同時(shí)C5Frac還可以避免近似解在濾失較大時(shí)出現(xiàn)的數(shù)值波動(dòng)。這表明在不同的濾失系數(shù)范圍下,所建立的方法均能準(zhǔn)確捕捉裂縫尺寸和濾失程度。

圖2 C5Frac解和參考解的裂縫半徑、縫口寬度和進(jìn)液效率對(duì)比(單裂縫擴(kuò)展驗(yàn)證)Fig.2 Comparison between the C5Frac solution and the reference solution in terms of fracture radius, width at fracture mouth and fluid efficiency (validation of single fracture propagation)

2.2 多裂縫擴(kuò)展

在證明單裂縫擴(kuò)展結(jié)果相對(duì)單裂縫參考解的準(zhǔn)確性后,進(jìn)一步驗(yàn)證C5Frac在模擬多裂縫擴(kuò)展方面的準(zhǔn)確性。選取ILSA II模型進(jìn)行對(duì)比驗(yàn)證,它是ILSA[18]的擴(kuò)展版本,一種基于隱式水平集算法的平面三維裂縫擴(kuò)展全耦合模型[17]。ILSA II使用三維位移不連續(xù)法求解彈性方程,使用有限體積法求解流體流動(dòng)。ILSA II的新穎之處在于可考慮裂縫尖端的多種漸進(jìn)行為,使用隱式水平集方法追蹤裂縫移動(dòng)邊界,即使在較粗糙的網(wǎng)格下仍能夠獲得極高的計(jì)算精度,適合作為驗(yàn)證多裂縫擴(kuò)展的參考解。

建立的模型共包含5簇,簇間距為30 m。模型相關(guān)參數(shù)設(shè)定如表1所示。

表1 模型參數(shù)設(shè)置

圖3為C5Frac和ILSA II模擬得到的多裂縫形態(tài)?？梢钥闯?在不同的模擬時(shí)間下,兩者的形態(tài)非常接近。同時(shí),還觀察到明顯的應(yīng)力陰影現(xiàn)象,外側(cè)裂縫擴(kuò)展具有優(yōu)勢(shì),而內(nèi)側(cè)和中間裂縫的擴(kuò)展受到抑制。

圖3 C5Frac和ILSA II在不同時(shí)間下的多裂縫形態(tài)Fig.3 Geometry of multiple fractures at different time between C5Frac and ILSA II

在此基礎(chǔ)上,進(jìn)一步統(tǒng)計(jì)各條裂縫的相關(guān)參數(shù)進(jìn)行定量對(duì)比,包括:裂縫無(wú)因次半徑(裂縫半徑與段長(zhǎng)的比值),縫內(nèi)流量和縫口寬度。統(tǒng)計(jì)結(jié)果顯示,兩者計(jì)算得到的無(wú)因次半徑、縫內(nèi)流量和縫口寬度非常接近。

從圖4(b)可以看出,在模擬15 s后,泵注排量幾乎全部分配給外側(cè)裂縫,流向內(nèi)側(cè)和中間裂縫的排量降至接近零。這是因?yàn)閮?nèi)側(cè)和中間裂縫在外側(cè)裂縫產(chǎn)生的誘導(dǎo)應(yīng)力場(chǎng)中相互競(jìng)爭(zhēng)擴(kuò)展,而內(nèi)側(cè)和中間裂縫的任何額外擴(kuò)展都會(huì)增強(qiáng)誘導(dǎo)應(yīng)力場(chǎng)。隨著模擬時(shí)間的增加,外側(cè)裂縫的擴(kuò)展變得更加明顯,而內(nèi)側(cè)和中間的尺寸幾乎不再發(fā)生變化。上述結(jié)果表明,本文方法能夠準(zhǔn)確模擬多裂縫競(jìng)爭(zhēng)擴(kuò)展過(guò)程,且具有極高的計(jì)算精度。

圖4 C5Frac解和參考解的裂縫無(wú)因次半徑、縫內(nèi)流量、縫口寬度和裂縫總面積對(duì)比(多裂縫擴(kuò)展驗(yàn)證)Fig.4. Comparison between the C5Frac solution and the reference solution in terms of dimensionless fracture radius, flow rate inside fracture, width at fracture mouth and total fracture area (validation of multiple fractures propagation)

由于壓裂后單井的產(chǎn)能與儲(chǔ)層改造程度高度相關(guān),因此在統(tǒng)計(jì)各條裂縫參數(shù)的基礎(chǔ)上,進(jìn)一步將裂縫總面積作為評(píng)價(jià)水力壓裂效果的重要指標(biāo),其定義為

(41)

在模擬的初期,由于所有裂縫的尺寸都較小,因此它們之間的互作用力不顯著,所有的裂縫都以較為接近的擴(kuò)展速率增大裂縫面積,并且該過(guò)程幾乎與時(shí)間呈線性關(guān)系。從圖4(d)中的對(duì)比曲線可看出,ILSA II求出的總裂縫面積隨時(shí)間變化的曲線與參考解的結(jié)果誤差保持在5%以內(nèi),這進(jìn)一步說(shuō)明C5Frac模擬的準(zhǔn)確性。

2.3 計(jì)算效率對(duì)比

隨著模擬對(duì)象尺度的增加,各種需要?jiǎng)澐志W(wǎng)格的方法(離散元,有限元,擴(kuò)展有限元)的計(jì)算時(shí)間都會(huì)大幅度增加。而本文模型不會(huì)隨著模擬時(shí)間的增加降低計(jì)算速度,其在保證精度的同時(shí)所擁有的計(jì)算效率優(yōu)勢(shì)是現(xiàn)有方法中獨(dú)一無(wú)二的。如表2所示,將C5Frac和ILSA II的計(jì)算時(shí)間進(jìn)行對(duì)比,可以發(fā)現(xiàn)當(dāng)模擬時(shí)間相同時(shí),C5Frac的計(jì)算時(shí)間卻遠(yuǎn)低于ILSA II,計(jì)算效率得到大幅提高。

表2 計(jì)算效率對(duì)比

3 結(jié)論

(1)將無(wú)量綱近似解與能量守恒相結(jié)合,建立一種半解析半數(shù)值的多裂縫擴(kuò)展模擬新方法,并開(kāi)發(fā)計(jì)算程序C5Frac。

(2)在進(jìn)行單裂縫擴(kuò)展模擬時(shí),C5Frac在裂縫半徑、縫口寬度和進(jìn)液效率方面的計(jì)算結(jié)果與參考解的誤差均在1%以內(nèi)。

(3)在進(jìn)行多裂縫擴(kuò)展模擬時(shí),C5Frac計(jì)算得到的裂縫半徑、縫口寬度和縫內(nèi)流量與ILSA II模型非常接近,裂縫總面積與ILSA II模型的誤差保持在5%以內(nèi)。

(4)在能夠保證計(jì)算準(zhǔn)確度的前提下,C5Frac的計(jì)算效率要遠(yuǎn)高于ILSA II,大幅減少模擬所需時(shí)間,能夠更加快速高效地模擬多裂縫擴(kuò)展。