考慮單元劈裂的流-固耦合連續(xù)-非連續(xù)方法及定向水力壓裂模擬

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

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

1 引 言

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

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

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

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

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

2 可沿對角線開裂的流-固耦合方法

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

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

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

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

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

2.2 流體對固體作用力的計算

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

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

(1)

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

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

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

(2)

fs=s2(3-2s)

(0≤s≤1) (3)

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

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

(4a)

(4b，4c)

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

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

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

(5)

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

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

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

(6a)

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

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

(6b)

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

2.3 計算流程

3 方法驗證

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

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

參數(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。時間步長Δt取為 5×10-7s。計算在平面應(yīng)變和大變形條件下進行，不計重力。為了與理論解對比，在計算時，裂縫寬度保持不變。這樣，流體對固體的作用力將不受固體參數(shù)的影響。

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

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

3.2 KGD模型

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

圖4 射孔尖端處裂縫張開度和裂縫半長隨時間的變化Fig.4 Variation of the fracture opening displacement at the perforation tip and half fracture length with time

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

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

表1 KGD模型的部分計算參數(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。為了便于對比，模型均剖分為5695個四邊形單元。在模型中預(yù)置射孔，長度約為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，其他計算參數(shù)取值列入表2。計算在平面應(yīng)變和大變形條件下進行，不計重力。

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

4.2 計算步驟和方案

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

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

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

4.3 射孔角度的影響

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

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

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

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

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

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

4.4 σx的影響

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

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

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

圖8 不同β時監(jiān)測位置處流體壓力及裂縫區(qū)段數(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時裂縫的時空分布Fig.9 Spatiotemporal distribution of fractures under different σx

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

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

4.5 數(shù)值與實驗結(jié)果對比

實驗?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時監(jiān)測位置處流體壓力及裂縫區(qū)段數(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給出了裂縫擴展路徑的實驗[2]和數(shù)值結(jié)果。可以看出，兩種結(jié)果在定性上吻合。由于本文的模擬條件為二維，不能考慮σz的影響，所以，裂縫轉(zhuǎn)向的數(shù)值結(jié)果與實驗結(jié)果尚存在一些差別。

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

5 結(jié) 論

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

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

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