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

白雪元， 王學濱， 劉桐辛

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

1 引 言

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

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

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

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

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

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

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

2.1 連通裂縫網絡的構建

欲計算流體在裂縫中的流動，需建立連通裂縫網絡。具體構建過程[18]如下。

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

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

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

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

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

(1)

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

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

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

(2)

fs=s2(3-2s)

(0≤s≤1) (3)

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

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

(4a)

(4b，4c)

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

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

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

(5)

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

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

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

(6a)

與端點B對應的單元節(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的流體壓力，裂縫右端為無滲透邊界。

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

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

圖3給出了不同時間t時裂縫中流體壓力的理論解[14]和數值解?？梢钥闯觯粫r刻，自左向右流體壓力逐漸減?。浑S著t的增加，同一位置的流體壓力增加，最終將達到PL。數值解與理論解吻合較好，這在一定程度上驗證了方法的正確性。

3.2 KGD模型

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

圖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個正方形單元，四個端面均為法向約束。在模型左側中間位置預置了射孔，射孔長度為0.05 m，向射孔中注入流體的流量q0取為5×10-6m2/s。Δt取為 1×10-6s，其他計算參數取值列入表1。計算在平面應變和大變形條件下進行，不計重力。

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

表1 KGD模型的部分計算參數Tab.1 Some calculation parameters of the KGD model

4 射孔角度及水平應力的影響

4.1 模型和參數

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

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

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

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

4.2 計算步驟和方案

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

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

為了研究β和σx的影響。共設計了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的最大主應力及裂縫時空分布，應當指出，本文壓裂生成的裂縫均為拉裂縫，簡稱為裂縫；黑色線段表示裂縫(包括預置的裂縫和壓裂后生成的裂縫)；圓點表示預置射孔尖端?？梢钥闯觯芽p從預置射孔的尖端起裂，呈彎曲狀向模型內部擴展，擴展方向逐漸趨于最小主應力方向。在裂縫擴展過程中，最大主應力始終集中在裂縫尖端，裂縫附近會出現一些次生裂縫，但不影響主裂縫的擴展。

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

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

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

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

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

4.4 σx的影響

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

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

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

圖8 不同β時監(jiān)測位置處流體壓力及裂縫區(qū)段數目隨時步數目的變化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時裂縫中流體壓力分布(時步數目=100000)Fig.10 Distribution of the fluid pressure in fractures under different σx (timestep=100000)

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

4.5 數值與實驗結果對比

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

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

圖12給出了裂縫擴展路徑的實驗[2]和數值結果?？梢钥闯?，兩種結果在定性上吻合。由于本文的模擬條件為二維，不能考慮σz的影響，所以，裂縫轉向的數值結果與實驗結果尚存在一些差別。

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

5 結 論

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

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

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