基于PHF-LSM-MT法巖石材料細觀邊界建模與水力裂縫擴展研究

李 明, 趙 岐, 陳 昭, 梁 力

(東北大學 資源與土木學院，遼寧 沈陽 110819)

水力壓裂技術作為油氣增產手段已被廣泛應用于實際工程并取得顯著效果.目前關于水力壓裂擴展規律的研究手段可分為：物理試驗、理論研究和數值模擬.其中物理試驗是探究巖石水力裂縫擴展特征、分布形態以及影響起裂壓力因素的有效方法.例如Zhang等[1]采用三軸水力壓裂試驗研究天然裂縫的空間位置和傾角對水力裂縫擴展形態的影響.

在理論研究方面，已有的方法多基于斷裂力學的相關研究成果，如PK模型、KGD模型[2]和PKN模型[3]等.物理試驗雖然可以較真實地模擬原位裂縫擴展形態，但成本較高，理論求解雖然可以得到精確解，但均建立在嚴格的假設基礎上.

數值模擬方法具有低成本高收益的特點，可以彌補物理試驗和理論研究的不足.例如：內聚區模型[4]、彌散裂縫模型[5]、相場法[6]等.此外不同方法的結合亦可解決水力壓裂的相關問題，如Yan等[7]采用有限元法耦合離散元法實現了水力壓裂過程的模擬，并通過物理試驗證明壓裂裂隙與地應力的關系.

巖石的宏觀力學特性會受到細觀結構特征的影響[8]，因此有必要考慮天然巖石細觀邊界特征對力學特性和水力壓裂過程的影響.如Li等[9]結合水平集法(level set method,LSM)建立了含有圓形、橢圓形以及空間非均質的巖石材料模型，并結合基于滲透率修正的水力壓裂模型(permeability-based hydraulic fracture, PHF)分析了水力裂縫的傳播過程.

已有的PHF-LSM法可以考慮顆粒分布對水力壓裂過程的影響，但其中未考慮積分點處細觀邊界特征.本文結合柏林噪聲和橢圓水平集方程，建立以顆粒半徑為參數的非規則顆粒建模方法，基于重心坐標法[10]實現非規則顆粒填充，并引入Mori-Tanaka(MT)均勻化方法[11]考慮高斯積分點處的細觀邊界特征.最終，建立PHF-LSM-MT法的非均質巖石材料的水力壓裂數值計算模型，并分析考慮細觀特征的非規則巖石材料的水力裂縫特性.

1 非規則顆粒邊界方程

在本文的研究中假設非規則顆粒Ω由顆粒包含的橢圓形(Ω0)和非規則圓環(Ω1)組成，如圖1所示.

圖1 非規則顆粒等效分解

極坐標系下的非規則顆粒邊界的水平集方程，如式(1)所示.

φi(θ)=r0(θ)+r1(θ) (0<θ≤2π).

(1)

2 顆粒填充和碰撞檢測

假設巖石內部的包裹體之間不重疊，在二維區域內生成多個非規則顆粒時，顆粒之間的關系可分為：包含、相交、相切和相離，如圖2所示.本文通過包裹體邊界上的點是否在其他包裹體中的方法來判別顆粒間的相對位置關系.

圖2 非規則顆粒位置關系

2.1 非規則顆粒接觸判定

將顆粒離散為三角形單元集合，邊界由線性差值簡化，如圖3a所示，其面積由式(2)計算.

圖3 非規則顆粒碰撞檢測方法

(2)

式中：Sj為第j個三角形的面積；N為三角形個數，本文建模過程中取N=360.

重心法原理如圖3b所示，直角坐標系下三角形O1O2O3的頂點坐標向量分別為O1，O2，O3，點P′在三角形內部的權重公式為[13]

P′-O1=u(O3-O1)+v(O2-O1).

(3)

式中：u，v均為常數，當u=0，v=0時,P′代表O1點；當u=0,v=1時,O′代表O2點；當u=1,v=0時,P′代表O3點.根據圖中P′，P1，P2，P3可知，如果判斷點在三角形內部，則必須要滿足①u≥0；②v≥0；③u+v≤1，取等號時，表示判斷點在三角形邊上.令V1=(P′-O1)，V2=(O3-O1)，V3=(O2-O1)，將式(3)的等式兩邊分別乘以V2和V3即可得到v和u.

2.2 非規則顆粒間的位置關系

建立非規則顆粒I1和I2并分別包含以下信息：中心坐標(x0,y0)；長半軸(a)；短半軸(b)；顆粒傾角(β)；圓環厚度均值(A)；顆粒邊界頂點坐標(xj,yj).xj和yj分別代表第j個頂點的x坐標和y坐標，其中1≤j≤360.頂點坐標的儲存順序按照逆時針方向，當長半軸與水平夾角為0°時記為第一個頂點坐標.

在調用顆粒信息時采用以下書寫格式(調用第1個顆粒的中心坐標記為I1(x0,y0))：

①顆粒I1包含顆粒I2：I2(x0,y0)，…，I2(x360,y360)均在I1內部，判別式恒成立.②顆粒I2包含顆粒I1：I1(x0,y0)，…，I1(x360,y360)均在I2內部，判別式恒成立.③顆粒I2與顆粒I1相交：I2(xk,yk)，…，I2(xl,yl)均在I1內部，共(l-k+1)個判別式成立；或I1(xm,ym)，…，I1(xn,yn)均在I2內部，共(n-m+1)個判別式成立，1≤k

當顆粒在進行碰撞判定時會出現圖4所示的4種特殊情況.通過細化顆粒邊界和控制最小距離Dmin避免圖4所示現象，其中Dmin與長半軸(a)和厚度均值(A)之間的關系見式(4).

圖4 非規則顆粒間特殊位置

Dmin+Dgap≥Iia+IiA(Dgap≥0).

(4)

式中：Dgap為顆粒間縫隙；Iia和IiA分別為新顆粒Ii的長半軸和厚度均值.

2.3 非規則顆粒間距計算

以圖5所示五邊形為例，在確定I2的中心點V0不在顆粒I1內部的基礎上，計算向量OV0與水平方向夾角θt，在顆粒I1中一定可以找到2個連續的頂點坐標Vj+1(xj+1,yj+1)和Vj(xj,yj)，該頂點坐標和I1中心點連線與水平方向的夾角分別為θj+1，θj，三者之間的關系滿足θj+1≥θt≥θj.OVt的長度采用式(5)和式(6)計算，其中Vt為OV0與Vj+1Vj的交點.

圖5 非規則顆粒間距計算

|OVj×OVj+1|=|OVj×OVt|+|OVt×OVj+1|,

(5)

(6)

Dmin=min{|OV0|-|OVt|,|V0V1|,…,|V0V360|}.

(7)

通過式(4)和式(7)可計算不同顆粒間的最小距離Dgap.采用Processing軟件在給定區域內填充不同幾何特征的非規則顆粒模型，如圖6所示，其中a=(2.0～30.0)mm，b=(0.0～1.0)×a，β=(0～2π).

圖6 非規則顆粒填充

3 多尺度非規則邊界有限元模型

已有的基于水平集法建立的有限元模型如圖7所示.由于積分點的有效區域內會被賦值為單一材料屬性，顆粒邊界會呈現“鋸齒狀”,如圖7b所示.

圖7 隱式建模方法

3.1 非均質巖石材料邊界水平集方程及參數均勻化

以圖8所示含有隨機分布非規則顆粒(Dgap>0)，邊長為l的巖石模型為例.該模型受到初始水壓力(pw)、水平和垂直方向應力(σh，σv)作用，模型底邊受到垂直方向約束，底邊中心位置為固定約束.模型中間施加注水荷載.由式(1)建立的水平集方程為

圖8 非均質巖石力學模型

I={x∈Ω:φi(x)≤0} (i=1,2,…,n),

(8)

φi=(x(r,θ))=φi(θ)-r.

(9)

式中：x為平面Ω內積分點的位置向量；當φi<0時表示積分點在顆粒內部，反之，當φi>0時表示積分點在顆粒外部，當φi=0時表示積分點在顆粒邊界上.

當積分點的有效區域wsub×hsub內含有材料邊界時(如圖7中積分點i1，i2，i5，i9)，已有的簡化方法是通過階躍函數、線性函數和非線性函數對材料進行過度處理[14]，該三種方法均不能考慮材料邊界特征的影響.本文采用均勻化方法對高斯積分點有效區域進行均勻化處理，以消除邊界“鋸齒狀”現象.均勻化后的等效材料參數Λ由式(10)表示，即

Λ=f(I0,I1,…,In).

(10)

其中：I0代表基質的材料參數；I1，…，In代表顆粒的材料參數.PHF模型中共考慮了4種材料參數，其中彈性參數(E,v)采用MT[15]法計算，即：

(11)

(12)

Tr={I+Sr:(Ir-I0)}-1(r=1,…,n).

(13)

式中：I和I0是四階單位張量和基質的有效剛度張量；Sr代表第r相顆粒的四階Eshelby張量[16].

文獻[17]的研究表明等效材料彈性參數均在Voigt上限(VR+)和Reuss下限(VR-)之間，即

(14)

(15)

式中：Ma=κ1/κ0；κ0和κ1分別為基質和顆粒的初始滲透系數；φ1為顆粒的體積分數；ζ是與顆粒形狀有關的參數，本文采用的非規則顆粒是以橢圓形為基礎，因此計算中ζ取0.5[19].

3.2 顆粒邊界細觀模型驗證

本文依托ABAQUS軟件進行數值模擬，通過用戶子程序實現圖9所示積分點有效區域參數均勻化.步驟如下：①通過子程序UEXTERNALDB讀取非規則顆粒參數信息，根據式(11)建立隱式材料邊界模型，見圖9a.②通過子程序SDVINI調用單元信息，根據積分點類型獲得坐標信息和有效區域范圍wsub×hsub，見圖9b.③將有效區域離散成wmin×hmin子網格并計算每個子網格中基質和顆粒占有效面積的體積分數φ0和φ1，見圖9c.本文取值wmin=wsub/10，hsub=hmin/10.④根據式(11)、式(14)和式(15)對材料參數均勻化，并重新賦值到積分點，見圖9d.⑤重復循環完成所有單元檢測.

圖9 積分點參數均勻化

通過建立顆粒邊界并對不同區域賦值材料屬性得到模型(顯式模型)，其結果作為對照組，驗證LSM-MT方法的可行性.以模型垂直方向楊氏模量Ey為例，對比LSM法與LSM-MT法建立模型(隱式模型、改進模型)相對于顯式模型計算結果的相對誤差er=|(Ξim-Ξex)/Ξex|和絕對誤差εr=Ξim-Ξex，其中Ξim為隱式方法計算結果，Ξex為顯式方法計算結果，其相對誤差和絕對誤差記為eri和εri，改進模型的相對誤差和絕對誤差記為erm和εrm.

三種方法建立的含有非規則顆粒分布的有限元模型如圖10所示，可以看出隱式模型中材料邊界呈“鋸齒狀”，而改進模型中材料邊界有明顯改善.同時采用表1所示包裹體的幾何參數，建立基質和顆粒的楊氏模量分別取值20 GPa和80 GPa，泊松比分別取值0.25和0.2的模型.通過不同方法計算的等效模量如圖11所示，結果表明，等效彈性參數Ey均落在VR上下界限內，且在MT法的解上下波動.三種模型的Ey均隨顆粒體積分數增加而增大.

圖10 基于不同建模方法有限元模型

表1 非規則顆粒幾何參數取值

圖11 體積分數改變對彈性模量的影響

兩種隱式模型在不同體積分數下的相對誤差和絕對誤差如圖12所示.

圖12 不同方法的有效模量誤差對比

結果表明兩種模型計算得到的彈性模量均大于理論值，且LSM-MT法的計算誤差均小于LSM法.因此本文建立的LSM-MT建模方法可有效考慮材料邊界分布特征，并增加計算結果的精度.

4 非均質巖石的水力裂縫特性

采用基于Biot固結理論建立的PHF水力壓裂數值計算模型，結合LSM-MT法建立含有非規則顆粒的巖石模型實現其水力壓裂過程的模擬.

4.1 PHF方法的控制方程

在深部巖層中水力壓裂過程是滲流與裂隙巖體相互作用的過程，其平衡方程、滲流方程、幾何方程、物理方程和有效應力原理如下[9]：

σij,j+ρfi=0 (i,j=1,2,3),

(16)

(17)

(18)

(19)

(20)

(21)

式中：參數κave=(κmax+κ0)/2，κmax和κ0分別為滲透系數上限值和初始值；κδ=(κmax-κ0)/2；Rm和σm分別為巖石材料的抗拉強度和平均有效應力；ξ為損傷系數；μmix為混合黏滯系數；krw為水相對滲透率，采用Corey[20]關系計算，即

krw=krw0((Sw-Swi)/(1-Swi-Sorw))Nw.

(22)

式中：krw0，Sw，Swi，Sorw和Nw分別為水相對滲透率端點、飽和度、束縛水飽和度以及剩余油飽和度；本文計算中krw0=0.5，Nw=3.66，Swi=Sorw=0.15[20].

4.2 PHF-LSM-MT水力壓裂數值計算模型

PHF-LSM-MT水力壓裂數值模型的建立采用ABAQUS軟件二次開發功能，流程如圖13所示.子程序SDVINI包含8個狀態變量：彈性模量(E)、泊松比(v)、基巖抗拉強度(Rm)、有效滲透系數(κ)、基巖密度(ρ)、飽和度(S)、孔隙比(e)和混合黏滯系數(μmix).依據式(11)實現非規則顆粒材料分布，并采用均勻化方法修改材料邊界處積分點有效參數.子程序USDFLD用來接受SDVINI傳入的狀態變量，并通過式(21)修正等效開裂區域滲透系數.

圖13 PHF-LSM-MT模型算法實現

4.3 均質巖石材料水力裂縫擴展

以邊長為100 mm的二維均質巖石模型為例，將PHF模型計算得到的水力裂縫擴展形態與KGD模型的理論解進行對比.力學模型如圖8所示，材料參數采用某地下3 088 m的實測數據，如表2所示.

表2 基巖材料力學參數

模型受到初始水平和垂直方向地應力為σh= 61.47 MPa和σv= 71.62 MPa，初始水壓力pw= 52 MPa，注水速度qinj=0.001 mm/s，注水時間為600 s.由于PHF模型是通過修改破壞區滲透系數表示裂縫擴展過程的一種彌散裂縫模型，因此模擬過程中沒有出現真實的裂縫，計算得到的等效開裂區域發展過程如圖14所示.采用PHF法計算的裂縫等效高度hf為開裂區域的垂直方向范圍，裂縫等效寬度wf為注水點位置裂縫開裂區域的水平方向位移差(ur-ul).

圖14 裂縫等效開裂發展過程

基于KGD理論模型的水力裂縫高度hf、裂縫開度wf和注水壓力pw的理論解如式(23)～式(25)[9]所示：

(23)

(24)

(25)

式中：G為剪切模量；當水力裂縫為雙向擴展時，參數a，b，c分別為0.48，1.32和1.91[9].

通過PHF模型計算的水力裂縫等效高度、等效開度和注水點水壓力隨注水時間變化過程與KGD理論解對比如圖15所示.

圖15 水力裂縫特征驗證

該計算結果表明PHF模型計算的數值解與解析解基本吻合.由于KGD模型假設裂縫處于無限半空間中，而在注水過程中隨著裂縫不斷擴展并逐漸接近模型邊界，導致兩者有一定偏差.

4.4 非均質巖石水力裂縫發展特性

建立圖8所示的非均質巖石的數值計算模型，邊界條件與均質巖石情況相同.計算中假設包裹體不易開裂，所用材料參數如表3所示.每種包裹體體積分數包括3個隨機算例.

表3 包裹體材料力學參數

當注水時間為600 s時，基于PHF-LSM-MT法計算的非均質巖石模型等效開裂區域如圖16所示.該結果表明：當傳播路徑存在包裹體時，水力裂縫在材料交匯處發生變化.當包裹體與裂縫發展方向的有效面積和傾角均較大，裂縫沿包裹體表面擴展，如圖16a和圖16c所示；當包裹體在裂縫發展方向的有效面積較大但傾角較小時，包裹體會阻止裂縫擴展，隨注水時間增加，裂縫沿包裹體表面發展，最終繞過包裹體，如圖16c和圖16d所示；當包裹體在裂縫發展方向的有效面積較小時，裂縫直接繞過包裹體繼續發展，如圖16b和圖16d所示.

圖16 非均質巖石水力裂縫模型等效開裂區域

注水點位置的水壓力變化曲線如圖17所示，兩種建模方法得到的結果均呈現先增加至起裂壓力后降低至傳播壓力的趨勢，與已有的計算結果一致[14].

圖17 注水點單元水壓力發展歷程

注水位置應力路徑變化如圖18所示.結果表明包裹體分布影響了應力路徑初始應力狀態.當注水時間增加時，平均有效應力p由初始狀態降低至零(受壓狀態)，隨后開始增加(受拉狀態).剪應力q在開裂過程處于降低狀態，最終為零.

圖18 注水點單元應力路徑發展歷程

圖19給出了不同體積分數，兩種方法計算的起裂壓力.當包裹體體積分數相同時，三種隨機工況均表明改進模型計算得到的起裂壓力高于隱式模型，絕對誤差最大值為0.101 MPa.基于線彈性理論的起裂壓力理論解為pb=3σh-σv+Rm= 116.962 MPa[9]，與本文中隱式模型和改進模型計算結果的相對誤差分別在8.98%～12.02%和8.94%～11.98%之間.

圖19 不同方法計算的起裂壓力

5 結 論

1)建立了基于柏林噪聲的非規則幾何邊界的水平集方程，采用重心坐標法檢測顆粒碰撞，實現以粒徑、長短半軸比、傾角和填充率為參數的非規則顆粒填充算法.

2)在已有LSM建模方法的基礎上，采用均勻化方法建立考慮非規則顆粒細觀邊界特征的LSM-MT方法.通過與顯式方法的計算結果對比，相對誤差為0.03%～2.58%，與LSM模型計算結果相比，改進幅度為0.12%～2.30%，結果表明LSM-MT法可以有效考慮非規則顆粒細觀邊界特征.

3)所建立的PHF-LSM-MT方法實現了非均質巖石材料建模及水力裂縫傳播過程模擬，通過與線彈性理論的起裂壓力對比，得到隱式模型和改進模型的相對誤差分別在8.98%～12.02%和8.94%～11.98%之間，且隨包裹體含量的增加，模型起裂壓力降低.