白雪元， 王學(xué)濱， 劉桐辛

(1.遼寧工程技術(shù)大學(xué) 力學(xué)與工程學(xué)院,阜新 123000；2.遼寧工程技術(shù)大學(xué) 計(jì)算力學(xué)研究所,阜新 123000)

1 引 言

水力壓裂過(guò)程中，裂縫的擴(kuò)展受到巖性、射孔角度和地應(yīng)力等影響[1]。在試驗(yàn)方面，文獻(xiàn)[2，3]采用相似材料模型研究了射孔角度對(duì)起裂壓力和裂縫擴(kuò)展路徑的影響。

目前，水力壓裂過(guò)程的數(shù)值模擬多是基于連續(xù)方法進(jìn)行的。其中，擴(kuò)展有限元方法適于模擬裂縫擴(kuò)展，且裂縫擴(kuò)展不受網(wǎng)格影響。盛茂等[4]采用擴(kuò)展有限元方法模擬了恒定水壓作用下水力裂縫的擴(kuò)展。趙金洲等[5]采用擴(kuò)展有限元方法研究了射孔角度對(duì)裂縫擴(kuò)展路徑和裂縫面剪應(yīng)力、切向位移的影響。連續(xù)方法相對(duì)成熟，已得到廣泛應(yīng)用，但難以處理多條任意分布的裂縫相互作用和大尺度材料的運(yùn)動(dòng)問(wèn)題[6]。

在非連續(xù)方法中，Zhang等[7]采用顆粒離散元方法研究了射孔角度對(duì)裂縫擴(kuò)展的影響。在顆粒離散元方法中，顆粒之間存在一定的孔隙。當(dāng)裂縫擴(kuò)展時(shí)，裂縫與孔隙混在一起，裂縫擴(kuò)展和流體流動(dòng)路徑并不直觀。Jiao等[8]基于非連續(xù)變形分析方法構(gòu)建了水-力耦合模型，模擬了不同射孔角度時(shí)巖樣的破裂過(guò)程。在非連續(xù)方法中，盡管可以通過(guò)在單元之間加入黏結(jié)單元來(lái)模擬破裂[9]，但通常需引入法向和切向接觸剛度，這將對(duì)連續(xù)介質(zhì)的應(yīng)力和應(yīng)變計(jì)算精度產(chǎn)生影響。

為了彌補(bǔ)連續(xù)方法和非連續(xù)方法的不足，科技人員已發(fā)展了一些連續(xù)-非連續(xù)方法[10-13]。如有限元/離散元方法(FEM/DEM)、連續(xù)-非連續(xù)單元方法(CDEM)。這些方法適于模擬巖石變形及裂縫擴(kuò)展過(guò)程。在此基礎(chǔ)上，科技人員又發(fā)展了一些流-固耦合方法[14-16]。如嚴(yán)成增等[14]在FEM/DEM方法中加入了流-固耦合計(jì)算模塊，提出了一種FDEM-Flow方法。王理想等[15]提出了一種CDEM和中心型有限體積法相結(jié)合的流-固耦合方法。

為了更真實(shí)地模擬水力壓裂過(guò)程中的巖石變形、裂縫擴(kuò)展及流體流動(dòng)，首先，以可沿四邊形單元對(duì)角線開(kāi)裂的拉格朗日元與離散元耦合的連續(xù)-非連續(xù)方法[17]為基礎(chǔ)，引入流體計(jì)算模塊，發(fā)展了一種流-固耦合方法；然后，通過(guò)與單裂縫非穩(wěn)態(tài)滲流模型和KGD模型的理論解進(jìn)行對(duì)比，驗(yàn)證了方法的正確性；最后，研究了射孔角度及水平應(yīng)力對(duì)裂縫擴(kuò)展的影響。

2 可沿對(duì)角線開(kāi)裂的流-固耦合方法

在可沿四邊形單元對(duì)角線開(kāi)裂的拉格朗日元與離散元耦合方法的基礎(chǔ)上，考慮了流體的作用，發(fā)展了一種流-固耦合方法。在該方法中，裂縫可沿單元邊界和四邊形單元對(duì)角線擴(kuò)展；假定流體僅能在裂縫中流動(dòng)，固體是非滲透的；流體的流動(dòng)滿足立方定律。該方法適于模擬流體作用下裂縫的彎曲擴(kuò)展，不需要進(jìn)行網(wǎng)格重構(gòu)，而且，在不降低時(shí)間步長(zhǎng)的情況下，在一定程度上降低了網(wǎng)格依賴性。

2.1 連通裂縫網(wǎng)絡(luò)的構(gòu)建

欲計(jì)算流體在裂縫中的流動(dòng)，需建立連通裂縫網(wǎng)絡(luò)。具體構(gòu)建過(guò)程[18]如下。

記當(dāng)前非連通裂縫(包括新生裂縫與原有非連通裂縫)的集合為M，取M中每一個(gè)元素J進(jìn)行判斷，若J與連通裂縫相連，則J為連通裂縫，并刪除M中的J元素。

如此循環(huán)，直到M中不再有連通的裂縫。若M中有剩余裂縫，則為非連通裂縫。

2.2 流體對(duì)固體作用力的計(jì)算

假設(shè)裂縫、端點(diǎn)、單元和節(jié)點(diǎn)之間的關(guān)系如 圖1 所示，裂縫1和裂縫2為沿單元邊界開(kāi)裂的裂縫；裂縫3為沿對(duì)角線開(kāi)裂的裂縫。為了便于表述和理解，稱裂縫兩端的流體節(jié)點(diǎn)為端點(diǎn)，每個(gè)端點(diǎn)對(duì)應(yīng)兩個(gè)固體節(jié)點(diǎn)，如裂縫1包含端點(diǎn)1和端點(diǎn)2，端點(diǎn)1對(duì)應(yīng)節(jié)點(diǎn)4和節(jié)點(diǎn)10，端點(diǎn)2對(duì)應(yīng)節(jié)點(diǎn)5和節(jié)點(diǎn)11。對(duì)于裂縫1，端點(diǎn)2和端點(diǎn)1之間的壓力差為

Δp=p2-p1+ρg(z2-z1)

(1)

式中p1和p2分別為端點(diǎn)1和端點(diǎn)2的流體壓力；z1和z2分別為端點(diǎn)1和端點(diǎn)2在重力方向的坐標(biāo)，端點(diǎn)的坐標(biāo)為該點(diǎn)對(duì)應(yīng)的兩個(gè)固體節(jié)點(diǎn)坐標(biāo)的平均；ρ為流體的密度；g為重力加速度。

圖1 裂縫、端點(diǎn)、單元和節(jié)點(diǎn)之間的關(guān)系Fig.1 Relationship among cracks,endpoints,elements and nodes

根據(jù)立方定律，由端點(diǎn)1流入端點(diǎn)2的流量q為[19]

(2)

fs=s2(3-2s)

(0≤s≤1) (3)

式中s為端點(diǎn)的飽和度，k為裂縫的滲透系數(shù)，k=1/(12μ)，μ為流體動(dòng)力黏度;L為端點(diǎn)1和端點(diǎn)2之間的距離，即裂縫長(zhǎng)度；a為裂縫張開(kāi)度，取裂縫兩端張開(kāi)度的平均值；fs為端點(diǎn)處與s有關(guān)的系數(shù)，fs 1和fs 2用于防止兩端點(diǎn)之間無(wú)壓力差時(shí)因重力產(chǎn)生的不合理流量，下標(biāo)1和下標(biāo)2分別代表端點(diǎn)1和端點(diǎn)2。在計(jì)算滲流時(shí)，a需大于0，所以，需設(shè)定裂縫最小張開(kāi)度amin，當(dāng)a

由于端點(diǎn)2不僅與裂縫1相連，還與裂縫2和裂縫3相連。這樣，可求得端點(diǎn)2的總流量Q=q1+q2+q3。當(dāng)前時(shí)步端點(diǎn)2的流體壓力p為

(4a)

(4b，4c)

式中p0為上一時(shí)步的流體壓力，K為流體的體積模量，V和V0分別為當(dāng)前時(shí)步和上一時(shí)步端點(diǎn)的體積，其為與端點(diǎn)相連裂縫體積之和的一半。

圖2 流體對(duì)單元作用力的施加Fig.2 Application of the fluid pressure on elements

若式(4a)計(jì)算的流體壓力為負(fù)，則說(shuō)明端點(diǎn)的流體不足。此時(shí)，將流體壓力設(shè)為0，并更新端點(diǎn)的飽和度為

(5)

式中s0為上一時(shí)步端點(diǎn)的飽和度。當(dāng)s<1時(shí)，根據(jù)式(5)更新端點(diǎn)的飽和度。當(dāng)s>1時(shí)，s=1。

對(duì)于一個(gè)裂縫(圖2)，計(jì)算出兩端點(diǎn)的流體壓力后，可確定作用于單元邊界的線性分布力，再將單元邊界一半的線性分布力等效到相應(yīng)單元節(jié)點(diǎn)上以獲得流體對(duì)相應(yīng)單元節(jié)點(diǎn)的作用力。若端點(diǎn)A和端點(diǎn)B的流體壓力分別為PA和PB，則與端點(diǎn)A對(duì)應(yīng)的單元節(jié)點(diǎn)3和單元節(jié)點(diǎn)5受到的流體作用力為

FA=(3PA/4+PB/4)×(L/2)

(6a)

與端點(diǎn)B對(duì)應(yīng)的單元節(jié)點(diǎn)4和單元節(jié)點(diǎn)6受到的流體作用力為

FB=(PA/4+3PB/4)×(L/2)

(6b)

若一個(gè)端點(diǎn)周圍有多個(gè)裂縫，則需將各裂縫中流體對(duì)固體節(jié)點(diǎn)的作用力累加以獲得流體對(duì)單元節(jié)點(diǎn)總的作用力。作用力的方向與裂縫角平分線方向(兩端點(diǎn)連線方向)垂直，并指向單元。

2.3 計(jì)算流程

3 方法驗(yàn)證

3.1 單裂縫非穩(wěn)態(tài)滲流模型

模型(圖3中插圖)長(zhǎng)W=1 m，寬H=0.1 m，剖分成100×10個(gè)正方形單元。模型中含有寬度為5×10-5m的裂縫。模型的上下兩個(gè)端面均為法向約束，在裂縫左端施加PL=9.5 MPa的流體壓力，裂縫右端為無(wú)滲透邊界。

參數(shù)包括固體參數(shù)和流體參數(shù)。固體一直處于彈性狀態(tài)，參數(shù)取值如下，面密度取為2400 kg/m2，彈性模量取為50 GPa，泊松比取為0.2，局部自適應(yīng)阻尼系數(shù)取為0.2。流體參數(shù)取值如下，μ=1×10-3Pa·s，K=20 MPa。時(shí)間步長(zhǎng)Δt取為 5×10-7s。計(jì)算在平面應(yīng)變和大變形條件下進(jìn)行，不計(jì)重力。為了與理論解對(duì)比，在計(jì)算時(shí)，裂縫寬度保持不變。這樣，流體對(duì)固體的作用力將不受固體參數(shù)的影響。

圖3 不同時(shí)刻裂縫中流體壓力分布Fig.3 Distribution of the fluid pressure in the fracture at different time

圖3給出了不同時(shí)間t時(shí)裂縫中流體壓力的理論解[14]和數(shù)值解?？梢钥闯?，同一時(shí)刻，自左向右流體壓力逐漸減小；隨著t的增加，同一位置的流體壓力增加，最終將達(dá)到PL。數(shù)值解與理論解吻合較好，這在一定程度上驗(yàn)證了方法的正確性。

3.2 KGD模型

KGD模型(圖4(a)中插圖)包括如下假定，(1) 裂縫高度固定；(2) 平面應(yīng)變； (3) 裂縫內(nèi)流體的流動(dòng)符合光滑平板流。

圖4 射孔尖端處裂縫張開(kāi)度和裂縫半長(zhǎng)隨時(shí)間的變化Fig.4 Variation of the fracture opening displacement at the perforation tip and half fracture length with time

模型尺寸為1 m×1 m，剖分成100×100個(gè)正方形單元，四個(gè)端面均為法向約束。在模型左側(cè)中間位置預(yù)置了射孔，射孔長(zhǎng)度為0.05 m，向射孔中注入流體的流量q0取為5×10-6m2/s。Δt取為 1×10-6s，其他計(jì)算參數(shù)取值列入表1。計(jì)算在平面應(yīng)變和大變形條件下進(jìn)行，不計(jì)重力。

圖4給出了射孔尖端處裂縫張開(kāi)度和裂縫半長(zhǎng)隨時(shí)間變化的數(shù)值解和理論解[15](圖4(b)中插圖為t=0.3 s和0.6 s時(shí)的垂直位移分布)。其中，裂縫半長(zhǎng)不包含射孔的長(zhǎng)度，從射孔尖端位置起裂開(kāi)始計(jì)時(shí)?？梢钥闯?，隨著t的增加，裂縫半長(zhǎng)增加；射孔尖端處裂縫張開(kāi)度呈先快后慢增長(zhǎng)。數(shù)值解與理論結(jié)果吻合較好，這進(jìn)一步驗(yàn)證了方法的正確性。

表1 KGD模型的部分計(jì)算參數(shù)Tab.1 Some calculation parameters of the KGD model

4 射孔角度及水平應(yīng)力的影響

4.1 模型和參數(shù)

定向射孔水力壓裂模型如圖5所示。模型尺寸為0.3 m×0.3 m，圓孔直徑為0.02 m。為了便于對(duì)比，模型均剖分為5695個(gè)四邊形單元。在模型中預(yù)置射孔，長(zhǎng)度約為0.03 m，射孔角度(射孔與x方向的夾角)為β。

圖5 定向射孔水力壓裂模型Fig.5 Hydraulic fracturing model containing directional perforations

在模型的左側(cè)面和上端面分別施加x方向水平應(yīng)力σx和y方向水平應(yīng)力σy=2 MPa，在模型的右側(cè)面和下端面施加法向約束。Δt取為1×10-7s，amin=2×10-5m，其他計(jì)算參數(shù)取值列入表2。計(jì)算在平面應(yīng)變和大變形條件下進(jìn)行，不計(jì)重力。

表2 水力壓裂模型部分計(jì)算參數(shù)Tab.2 Some calculation parameters of the hydraulic fracturing model

4.2 計(jì)算步驟和方案

(1) 施加載荷和約束后進(jìn)行計(jì)算，直至靜力平衡。此步驟共用4000個(gè)時(shí)步。

(2) 向射孔中注入流體，q0取為5×10-4m2/s，進(jìn)行計(jì)算。

為了研究β和σx的影響。共設(shè)計(jì)了7個(gè)計(jì)算方案。方案1～方案4的β分別為0°,30°,50°和70°，σx=2.5 MPa；方案5～方案7的σx分別為 1.5 MPa,2 MPa和3 MPa，β=50°。為了獲得破裂壓力，監(jiān)測(cè)了射孔靠近圓孔位置處(圖5)的流體壓力。

4.3 射孔角度的影響

圖6給出了方案3的最大主應(yīng)力及裂縫時(shí)空分布，應(yīng)當(dāng)指出，本文壓裂生成的裂縫均為拉裂縫，簡(jiǎn)稱為裂縫；黑色線段表示裂縫(包括預(yù)置的裂縫和壓裂后生成的裂縫)；圓點(diǎn)表示預(yù)置射孔尖端?？梢钥闯?，裂縫從預(yù)置射孔的尖端起裂，呈彎曲狀向模型內(nèi)部擴(kuò)展，擴(kuò)展方向逐漸趨于最小主應(yīng)力方向。在裂縫擴(kuò)展過(guò)程中，最大主應(yīng)力始終集中在裂縫尖端，裂縫附近會(huì)出現(xiàn)一些次生裂縫，但不影響主裂縫的擴(kuò)展。

圖6 方案3的最大主應(yīng)力及裂縫時(shí)空分布(單位：Pa)Fig.6 Spatiotemporal distribution of the maximum principal stress and fractures in scheme 3 (unit:Pa)

圖7給出了不同β時(shí)流體壓力的時(shí)空分布(數(shù)字表示方案編號(hào))?？梢钥闯?，當(dāng)β不同時(shí)，裂縫均從射孔尖端起裂；β對(duì)裂縫初始擴(kuò)展有明顯影響,β越大，裂縫發(fā)生轉(zhuǎn)向越明顯，轉(zhuǎn)向距離越大；裂縫最終總是偏向最小主應(yīng)力方向擴(kuò)展。這與文獻(xiàn)[2]通過(guò)實(shí)驗(yàn)獲得的β對(duì)裂縫擴(kuò)展的影響規(guī)律一致。在射孔附近，流體壓力最大；距離射孔越遠(yuǎn)，流體壓力越小。隨著時(shí)步數(shù)目的增加，裂縫中流體壓力降低，這與裂縫寬度增加有關(guān)。

圖8給出了不同β時(shí)監(jiān)測(cè)位置處流體壓力及裂縫區(qū)段數(shù)目隨時(shí)步數(shù)目的變化。應(yīng)當(dāng)指出，裂縫區(qū)段是指兩個(gè)單元之間的一段裂縫，多個(gè)裂縫區(qū)段依次相連形成裂縫。

從圖8(a)可以看出，隨著β的增加，射孔尖端起裂時(shí)，監(jiān)測(cè)位置的壓力呈增加趨勢(shì)；當(dāng)裂縫擴(kuò)展時(shí)，射孔位置處的壓力增加，這意味著壓裂難度增加。因此，在壓裂時(shí)，β與最小主應(yīng)力方向更接近有利于降低壓裂成本，提高壓裂效率，這與常識(shí)相符。隨著時(shí)步數(shù)目的增加，裂縫中流體壓力降低，這與圖7的流體壓力云圖結(jié)果一致。

從圖8(b)可以看出，隨著時(shí)步數(shù)目的增加，不同β時(shí)，裂縫區(qū)段數(shù)目的增速由快變慢，這說(shuō)明裂縫擴(kuò)展的速度逐漸變慢。這是由于隨著時(shí)間的增加，裂縫的體積增加變快，而流體注入的流量是恒定的，這使得裂縫中流體壓力逐漸降低，裂縫擴(kuò)展變慢。方案1～方案3的裂縫區(qū)段數(shù)目-時(shí)步數(shù)目曲線相近，而方案4的裂縫區(qū)段-時(shí)步數(shù)目曲線偏高。

4.4 σx的影響

圖9給出了不同σx時(shí)裂縫的時(shí)空分布?？梢钥闯觯瑑蓚€(gè)方向的水平應(yīng)力之差越大，裂縫偏轉(zhuǎn)距離越小。如方案7的兩個(gè)水平應(yīng)力之差大于方案3，而方案7的裂縫偏轉(zhuǎn)距離小于方案3。

圖10給出了不同σx時(shí)流體壓力分布(時(shí)步數(shù)目=100000)?？梢钥闯觯S著σx的增加，流體壓力的最大值增加。

圖7 不同β時(shí)裂縫中流體壓力時(shí)空分布(單位：Pa)Fig.7 Spatiotemporal distribution of the fluid pressure in fractures for different β (unit:Pa)

圖8 不同β時(shí)監(jiān)測(cè)位置處流體壓力及裂縫區(qū)段數(shù)目隨時(shí)步數(shù)目的變化Fig.8 Variation of the fluid pressure at the monitored position and the number of fracture segments with timesteps for different β

圖9 不同σx時(shí)裂縫的時(shí)空分布Fig.9 Spatiotemporal distribution of fractures under different σx

圖10 不同σx時(shí)裂縫中流體壓力分布(時(shí)步數(shù)目=100000)Fig.10 Distribution of the fluid pressure in fractures under different σx (timestep=100000)

圖11給出了不同σx時(shí)監(jiān)測(cè)位置處流體壓力及裂縫區(qū)段數(shù)目隨時(shí)步數(shù)目的變化規(guī)律?？梢钥闯觯S著σx的增加，起裂壓力增加，裂縫擴(kuò)展過(guò)程中監(jiān)測(cè)位置處的流體壓力增加(圖11(a))，也就是說(shuō)，裂縫擴(kuò)展時(shí)的壓裂難度增加。這與圖10的流體壓力云圖結(jié)果一致。隨著時(shí)步數(shù)目的增加，監(jiān)測(cè)位置處的流體壓力降低(圖11(a))。隨著時(shí)步數(shù)目的增加，裂縫的擴(kuò)展速度逐漸變慢(圖11(b))。σx對(duì)于裂縫區(qū)段數(shù)目演化規(guī)律的影響不明顯。

4.5 數(shù)值與實(shí)驗(yàn)結(jié)果對(duì)比

實(shí)驗(yàn)?zāi)Ｐ蜅l件[2]為σx=6 MPa，σy=1 MPa，σz=15 MPa，β=60°。應(yīng)當(dāng)指出，模型中含有套筒。據(jù)此建立二維模型，Δt取為5×10-8s，套筒面密度取為7800 kg/m2，彈性模量取為200 GPa，泊松比取為0.2，法向接觸剛度取為2×1013Pa/m，其他參數(shù)取值等均與4.1節(jié)相同。

圖11 不同σx時(shí)監(jiān)測(cè)位置處流體壓力及裂縫區(qū)段數(shù)目隨時(shí)步數(shù)目的變化Fig.11 Variation of the fluid pressure at the monitored position and the number of fracture segments with timesteps under different σx

圖12給出了裂縫擴(kuò)展路徑的實(shí)驗(yàn)[2]和數(shù)值結(jié)果?？梢钥闯觯瑑煞N結(jié)果在定性上吻合。由于本文的模擬條件為二維，不能考慮σz的影響，所以，裂縫轉(zhuǎn)向的數(shù)值結(jié)果與實(shí)驗(yàn)結(jié)果尚存在一些差別。

圖12 裂縫擴(kuò)展路徑的實(shí)驗(yàn)和數(shù)值結(jié)果Fig.12 Experimental and numerical results of fracture propagation paths

5 結(jié) 論

(1) 發(fā)展了一種單元劈裂的流-固耦合方法。在該方法中，裂縫可沿四邊形單元對(duì)角線和單元邊界擴(kuò)展，流體在裂縫中的流動(dòng)滿足立方定律。該方法適于模擬流體驅(qū)動(dòng)下的裂縫擴(kuò)展。

(2) 對(duì)于定向射孔水力壓裂模型，距離射孔越遠(yuǎn)，流體壓力越小。隨著時(shí)間的增加，裂縫中流體壓力降低。隨著射孔角度和x方向水平應(yīng)力的增加，裂縫起裂和擴(kuò)展過(guò)程中的流體壓力增加。兩個(gè)方向的水平應(yīng)力之差越大，裂縫轉(zhuǎn)向距離越小。

(3) 對(duì)于定向射孔水力壓裂模型，在裂縫擴(kuò)展過(guò)程中，裂縫區(qū)段數(shù)目的增速變慢，這與裂縫體積增加變快有關(guān)。