基于能量守恒的半解析半數值多裂縫擴展模擬方法

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

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

水平井分段多簇壓裂是開發油氣資源的有效增產措施之一,研究壓裂過程中水力裂縫的擴展延伸對油氣的高效開采是至關重要的[1],這直接關系到油氣井的產量。

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

為此,提出一種將無量綱近似解與能量守恒相結合的裂縫擴展模擬方法。該方法可充分考慮水力壓裂過程中的濾失影響及多裂縫間的競爭擴展效應。在此基礎上,編制半解析半數值裂縫擴展程序C5Frac,并將計算結果與ILSA II[17-18]的模擬結果進行對比驗證。該方法可以實現多裂縫擴展的快速模擬,且具有極高的計算精度,能夠為壓裂方案設計提供實時性指導。

1 半解析半數值裂縫擴展模型

1.1 控制方程

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

圖1 單個壓裂段內多裂縫擴展示意圖Fig.1 Schematic of multiple fractures propagation in a fracture stage

(1)

裂縫面的徑向擴展由井筒處注入的不可壓縮牛頓流體所驅動,并在具有滲透性的線彈性巖石中準靜態延伸(即遠低于巖石聲速),其特征為

(2)

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

韌性K′與巖石的斷裂韌性KIC的關系可表示為

(3)

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

提出上述假設之后,對于每條裂縫,需要求解的未知量包括:裂縫寬度w(r,t)、縫內壓力pf(r,t)、縫內體積流量q(r,t)、裂縫半徑R(t)、其他裂縫施加的作用力σ(r,t)以及縫口流量Q(t)。其中,r為裂縫內某一點到裂縫中心的距離,t為時間。控制方程如下。

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

(4)

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

(2)縫內流動方程由泊肅葉流動方程描述為

(5)

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

(6)

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

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

(7)

(4)根據線彈性斷裂力學,裂縫的Ⅰ型應力強度因子為[19]

(8)

當應力強度因子達到Ⅰ型斷裂韌性KIC時,裂縫即向前擴展。

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

(9)

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

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

(10)

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

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

整個系統的初始條件為

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)

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

(14)

此外,還能得到雷諾潤滑方程為

(15)

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

1.2 近似解

1.2.1 壓力分布近似解

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

(16)

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

(17)

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

1.2.2 韌性近似解

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

(18)

(19)

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

1.2.3 相互作用力近似解

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

(20)

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

相互作用力模型通過對所有相鄰裂縫的互作用應力求和來實現。因此,施加在水力裂縫上的互作用力近似表示為

(21)

1.2.4 彈性近似解

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

(22)

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

(23)

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

(24)

1.2.5 全局體積平衡

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

(25)

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

(26)

1.2.6 縫口條件

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

1.3 全局能量守恒

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

(27)

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

(28)

Df為與黏性流體流動相關的耗散率,可表示為

(29)

DL為與濾失相關的流體損失率,可表示為:

(30)

U為通過變形巖石應變能來增加應變能的一部分,是可恢復的彈性能,可表示為

(31)

(32)

Wo為原地應力對裂縫的作用。為考慮流體損失,原位應力的功被修正為

(33)

式(33)中:S為用于量化非均勻的原位應力;α為濾失相關系數。

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

(34)

通過解壓力、寬度和半徑對流量的隱性依賴表達式代替未知的μc和γ得

(35)

(36)

(37)

(38)

(39)

(40)

1.4 模型求解流程

求得近似解的步驟如下。

步驟1定義輸入參數σc、CL、E′、μ′、KI、Qo、h。

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

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

步驟6將時間步長向前推進Δt,Q(t)和α(t)可作為Q(t+Δt)和α(t+Δt)的預估值。

步驟7重復步驟3～步驟6,直至達到所需的總泵送時間。

2 模型驗證及計算效率對比

為驗證所建立的半解析半數值方法,將該方法模擬得到的計算結果與參考解結果進行對比。根據參考解所具有的適用性,采用兩類不同的參考解分別對單裂縫擴展和多裂縫擴展進行驗證,并對比多裂縫擴展模擬時的計算效率。

2.1 單裂縫擴展

單裂縫擴展的參考解是由Dontsov[24]開發的一種數值近似解計算得出的。該方法可捕捉水力裂縫在低濾失與高濾失間以及高黏度與強韌性狀態間的過渡行為,還可捕捉以韌性或黏度為主的水力裂縫和以濾失和存儲為主的水力裂縫的擴展行為,是一種可考慮濾失的裂縫擴展高精度模擬方法。

圖2為C5Frac模擬結果與參考解的對比結果,可以看出,兩者的誤差在1%以內,同時C5Frac還可以避免近似解在濾失較大時出現的數值波動。這表明在不同的濾失系數范圍下,所建立的方法均能準確捕捉裂縫尺寸和濾失程度。

圖2 C5Frac解和參考解的裂縫半徑、縫口寬度和進液效率對比(單裂縫擴展驗證)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 多裂縫擴展

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

建立的模型共包含5簇,簇間距為30 m。模型相關參數設定如表1所示。

表1 模型參數設置

圖3為C5Frac和ILSA II模擬得到的多裂縫形態。可以看出,在不同的模擬時間下,兩者的形態非常接近。同時,還觀察到明顯的應力陰影現象,外側裂縫擴展具有優勢,而內側和中間裂縫的擴展受到抑制。

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

在此基礎上,進一步統計各條裂縫的相關參數進行定量對比,包括:裂縫無因次半徑(裂縫半徑與段長的比值),縫內流量和縫口寬度。統計結果顯示,兩者計算得到的無因次半徑、縫內流量和縫口寬度非常接近。

從圖4(b)可以看出,在模擬15 s后,泵注排量幾乎全部分配給外側裂縫,流向內側和中間裂縫的排量降至接近零。這是因為內側和中間裂縫在外側裂縫產生的誘導應力場中相互競爭擴展,而內側和中間裂縫的任何額外擴展都會增強誘導應力場。隨著模擬時間的增加,外側裂縫的擴展變得更加明顯,而內側和中間的尺寸幾乎不再發生變化。上述結果表明,本文方法能夠準確模擬多裂縫競爭擴展過程,且具有極高的計算精度。

圖4 C5Frac解和參考解的裂縫無因次半徑、縫內流量、縫口寬度和裂縫總面積對比(多裂縫擴展驗證)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)

由于壓裂后單井的產能與儲層改造程度高度相關,因此在統計各條裂縫參數的基礎上,進一步將裂縫總面積作為評價水力壓裂效果的重要指標,其定義為

(41)

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

2.3 計算效率對比

隨著模擬對象尺度的增加,各種需要劃分網格的方法(離散元,有限元,擴展有限元)的計算時間都會大幅度增加。而本文模型不會隨著模擬時間的增加降低計算速度,其在保證精度的同時所擁有的計算效率優勢是現有方法中獨一無二的。如表2所示,將C5Frac和ILSA II的計算時間進行對比,可以發現當模擬時間相同時,C5Frac的計算時間卻遠低于ILSA II,計算效率得到大幅提高。

表2 計算效率對比

3 結論

(1)將無量綱近似解與能量守恒相結合,建立一種半解析半數值的多裂縫擴展模擬新方法,并開發計算程序C5Frac。

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

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

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