基于近場動力學方法的水力壓裂過程數(shù)值模擬

吳凡++李書卉++段慶林++李錫夔++張洪武

摘要： 利用近場動力學方法便于處理多裂紋萌生、擴展和分叉的優(yōu)點，將其應用于頁巖水力壓裂過程的數(shù)值模擬.結(jié)合頁巖水力壓裂機理提出在近場動力學中由破壞度跟蹤裂紋擴展路徑，并通過在新生成裂紋面法線方向施加水壓載荷的裂紋追蹤方法，成功模擬單射孔橫剖面開裂水壓實驗以及單射孔和多射孔水平井縱剖面的水力壓裂過程，得到壓裂導致的裂紋縫網(wǎng)結(jié)構(gòu).數(shù)值結(jié)果還表明初始射孔裂紋會顯著影響后續(xù)的水力壓裂過程.

關(guān)鍵詞： 頁巖氣開采； 水力壓裂； 近場動力學； 裂紋擴展； 多裂紋； 裂紋追蹤

中圖分類號： TE357.1文獻標志碼： B

Numerical simulation of hydraulic fracturing process

based on peridynamics method

WU Fan1， LI Shuhui1， DUAN Qinglin1， LI Xikui1， ZHANG Hongwu1

（State Key Laboratory of Structural Analysis for Industrial Equipment，

Dalian University of Technology， Dalian 116024， Liaoning， China）

Abstract： Due to the advantages of peridynamics method in dealing with the initiation， propagation and branching of multiple cracks， it is introduced into the numerical simulation on the hydraulic facture process of shale. Combining with the hydraulic fracture mechanism， a method for tracing crack path according to damage degree in peridynamics is proposed in which the hydraulic pressure is applied in the normal direction of the newly generated crack surface. The cracking pressure of single horizontal well in cross section and the hydraulic fracture of the horizontal well with single and multiple perforations in longitudinal section are successfully simulated by the method. The crack net due to fracturing is obtained. The numerical results also show that the initial perforations have significant effect in the subsequent hydraulic fracture process.

Key words： shale gas extraction； hydraulic fracture； peridynamics； crack propagation； multiple cracks； crack tracking

收稿日期： 2016[KG*9〗09[KG*9〗27修回日期： 2016[KG*9〗11[KG*9〗23

基金項目： 國家自然科學基金（11232003，11372066）；中央高?；究蒲袠I(yè)務費專項資金（DUT15LK07）；

遼寧省教育廳重點實驗室基礎(chǔ)研究項目（LZ2014002）

作者簡介： 吳凡（1992—），男，吉林遼源人，博士研究生，研究方向為近場動力學，（Email）251386704@qq.com

通信作者： 李書卉（1979—），女，遼寧沈陽人，博士，工程師，研究方向為材料破壞分析與實驗、并行計算技術(shù)等，（Email）shuhuili@dlut.edu.cn0引言

頁巖氣是蘊含在地下深處頁巖層區(qū)域有機質(zhì)泥頁巖內(nèi)的氣體，是以吸附和游離狀態(tài)為主的非常規(guī)天然氣資源.美國對頁巖氣的開采已有80多年的歷史.我國頁巖氣儲量位居世界第一，但開采起步相對較晚.[13]高效的頁巖氣開采技術(shù)對保證能源安全具有重要意義.國家“十二五”和“十三五”規(guī)劃均對頁巖氣開采產(chǎn)業(yè)進行強有力扶持.

由于頁巖埋藏深、質(zhì)地硬，頁巖氣被困于頁巖體內(nèi)難以逸出.目前，水力壓裂是最有效的頁巖氣開采技術(shù).該技術(shù)將壓裂液注入頁巖氣井中，通過高強度水壓壓裂巖石進而形成裂縫網(wǎng)絡，從而使得頁巖氣能通過由裂縫形成的高導流通道得以快速析出.

對水力壓裂過程實施準確有效的數(shù)值模擬有助于深入理解各工藝參數(shù)的作用機理，對裂縫形態(tài)進行預測，合理配置工藝參數(shù)，從而降低生產(chǎn)成本，提高經(jīng)濟效益.然而，由于該過程涉及多裂紋甚至海量裂紋的萌生、擴展和分叉，傳統(tǒng)模擬裂紋擴展的數(shù)值方法，如擴展有限元法等[45]，難以對其實施有效的數(shù)值模擬.

SILLING[6]提出近場動力學模型，其基本方程采用積分形式，因此在連續(xù)和不連續(xù)處均自然成立，受到力學界的廣泛關(guān)注[79].由于對裂紋無須進行特殊處理，近場動力學方法十分便于處理多裂紋尤其是海量裂紋的萌生、擴展和分叉問題.因此，該方法十分適合于頁巖水力壓裂過程的數(shù)值模擬，本文將對此進行初步探究.

1近場動力學理論

1.1連續(xù)模型

近場動力學基本方程基于牛頓第二定律的積分形式描述物質(zhì)點運動.假設結(jié)構(gòu)內(nèi)任一物質(zhì)點與其一定范圍內(nèi)的其他物質(zhì)點具有一定的相互作用，其基本方程為ρu··（x，t）=∫Hf（u（x^，t）-u（x，t），x^-x）dx^+

b（x，t） （1）式中：ρ，u和 b分別為密度、位移和外力體力項；f為力密度函數(shù)，被稱作對點力函數(shù)，代表物質(zhì)點x^對其臨近的物質(zhì)點x的作用；H為影響域，與無網(wǎng)格法的影響域定義相似，用于描述該物質(zhì)點x所能作用的其他物質(zhì)點所涵蓋的區(qū)域.

可視近場動力學為一種宏觀的分子動力學[10]，兩物質(zhì)點間的相互作用f可比擬為分子動力學中原子之間的范德華力，稱為鍵.對點力函數(shù)包含物質(zhì)本構(gòu)的所有信息，包括線性以及非線性等，因此通過修改對點力函數(shù)，對彈性、塑性、非線性問題均可以適用.對于微彈性材料，其對點力函數(shù)能夠由一個標量的微勢能函數(shù)的微分得到，即f（η，ξ）=ω（η，ξ）η （2）式中：ω（η，ξ）為微勢能函數(shù)，其物理意義為單個鍵內(nèi)儲存的能量，ξ=x^-x和η=u^-u分別為參考位置向量和相對位移向量.對于微彈性材料，有ω（η，ξ）=c（ξ）s2ξ2 （3）式中：ξ=ξ；η=η；s為鍵伸長，類似于應變，代表鍵的伸長比率.s定義為s=（ζ-ξ）/ξ （4） 式中：ζ=ξ+η為當前構(gòu)型下的相對位置向量，模量ζ=ξ+η代表當前鍵長.微模量c（ξ）表征材料剛度，物理意義為微觀下鍵的剛度.對于微彈性材料，c（ξ）與ξ無關(guān).在各項均勻拉伸變形下，將兩者應變能密度等效，可得出微模量與經(jīng)典力學中物理量的關(guān)系[6]為c=6Eπδ4（1-2ν） （5）式中：E為彈性模量；ν為泊松比；δ為影響域半徑.對于二維平面應力有c=9Eπδ3（1-2ν） （6）1.2斷裂模型

由動量守恒定律可知f平行于ζ，因此微彈性材料對點力函數(shù)可寫為f=ξ+ηξ+ηcsμ=ξ+ηξ+ηf （7）式中：f為對點力函數(shù)f的模量.因此，f=csμ （8）式中：μ為01函數(shù)，用以引入斷裂，表示該鍵的斷裂或完好，μ=1，s

0，其他 （9）式中：s0為極限伸長率，代表引發(fā)鍵斷裂的最大鍵伸長，當達到該設定值時，該鍵斷裂.微觀下鍵的斷裂在宏觀上能夠反映裂紋的發(fā)展情況.為在宏觀上描述斷裂，定義破壞度φ用以表示物質(zhì)點斷裂情況，φ（x，t）=1-∫H μ（x，t，ξ）dVH∫HdVH （10）式中：φ=0代表與該物質(zhì)點相連的所有鍵都保持完好，φ=1代表該點游離于結(jié)構(gòu).

2水力壓裂的近場動力學模型

水力壓裂過程示意見圖1.通過從地表鉆入豎井，穿過地下水層，深入到地下頁巖區(qū)域，然后換鉆頭打大斜度的水平井.水平井完成以后，通過送入井內(nèi)的射孔槍進行井內(nèi)射孔.射孔彈穿過井壁打入巖層，從而連通井內(nèi)和巖區(qū)形成通道.由于射孔彈穿透效果有限，需要在射孔之后進行水力壓裂.通過地面的壓裂車的高壓泵組將壓裂液注入水平井內(nèi)，利用高強度水壓壓裂頁巖儲層基質(zhì)形成復雜裂紋網(wǎng)絡，使裂縫壁面與儲層基質(zhì)的接觸面積更大，縮小油氣從任意方向的基質(zhì)向裂縫輸出的滲流距離，為頁巖氣形成高導流通道，加快頁巖氣的析出.[2]

Fig.1Schematic of hydraulic fracture process

水力壓裂是頁巖氣開采的主要增產(chǎn)措施，壓裂的體積改造是目前研究的熱點，即進一步提高儲層總剖面的使用程度，追求裂紋數(shù)目更多、裂紋形態(tài)更復雜、延伸距離更遠的縫網(wǎng)壓裂技術(shù).因此，本文針對水平井的橫縱剖面分別進行二維近場動力學數(shù)值模擬，對縫網(wǎng)的形成機制進行機理探究.

頁巖巖體物理參數(shù)取自吉林省農(nóng)安油田的頁巖實驗數(shù)據(jù)，見表1.近場動力學二維微彈性模擬的泊松比固定為1/3[10]，且由于泊松比對裂紋路徑、分叉的影響很小[11]，所以計算泊松比也取1/3.

Tab.1Parameters of shale materialE/GPaρ/（kg/m3）νs01.9032 4001/30.009 148 7

時間積分格式采用VelocityVerlet算法，射孔初始裂紋采用預先在結(jié)構(gòu)內(nèi)添設的方式實現(xiàn).為施加水壓，須對裂紋進行追蹤.為簡便起見，在本文的實現(xiàn)中，破壞度為0.35<φtrace<0.50的粒子被認為處于裂紋面上，而水壓的施加方向為斷裂鍵的“平均”方向.這里的“平均”是指對各斷裂鍵的取向根據(jù)其鍵長進行加權(quán)平均，其中權(quán)重系數(shù)αb（ξ）的計算式為αb（ξ）=∞i=1piξΔxi-1（11）式中：pi為擬合的多項式系數(shù).對于本文采用的均勻節(jié)點分布，在離散δ=3時權(quán)重系數(shù)αb（ξ）的取值見表2.

Tab.2Fitting points for directional weight coefficientsξ/Δx1223αb（ξ）1.001.000.300.01

3數(shù)值算例

通過幾個數(shù)值算例考察建立的近場動力學方法模擬水力壓裂過程的效果.為驗證近場動力學模型對裂紋擴展模擬的有效性和準確性，算例1模擬KalthoffWinkler實驗，算例2考察初始裂紋對所需壓裂水壓的影響情況，算例3考察單射孔水平井縱剖面的水力壓裂裂紋擴展情況，算例4進行多射孔水平井縱剖面的模擬.

3.1KalthoffWinkler實驗

KalthoffWinkler實驗是研究斷裂擴展角度的著名實驗.實驗中測得裂紋擴展角度呈70°.現(xiàn)采取的結(jié)構(gòu)尺寸見圖2.其中，v0=20 m/s.近場動力學模擬結(jié)果見圖3.

Fig.2KalthoffWinkler experiment setup

Fig.3KalthoffWinkler experiment simulation result obtained by peridynamics method

模擬得到的裂紋擴展角度與實驗測得的角度非常接近.裂紋在擴展一段距離之后便不再擴展（由于整體運動的產(chǎn)生），這一點與實驗結(jié)果也吻合.這些數(shù)值結(jié)果表明近場動力學方法能夠準確有效地模擬裂紋的擴展過程.

3.2初始裂紋對壓裂水壓的影響

該算例考慮水平井的橫剖面，分為有初始裂紋的套管井和無初始裂紋的裸眼井.水平井橫剖面示意圖和水壓加載見圖4.結(jié)構(gòu)尺寸為10.00 m×10.00 m，均勻離散下網(wǎng)格尺寸為Δx=0.02 m，節(jié)點數(shù)為65 000個，中央處水平井直徑為0.60 m，其中，套管井的射孔長度為0.60 m，裸眼井無射孔，對比2種井各自所需的起裂水壓.水壓加載見圖4b，增量為1 MPa/ms，0.200 ms后恒定（足夠大）.2種井的壓裂效果云圖見圖5.由此可看出：裸眼井的破裂壓力為13.2 MPa，而套管井的射孔在5.8 MPa時開始發(fā)生裂紋擴展，0.250 ms時發(fā)生分叉.該結(jié)果顯示射孔套管井能夠顯著降低所需的起裂水壓.

a）整體示意b）裸眼井與套管井示意c）水壓加載

Fig.4Schematic of cross section of horizontal

well and hydraulic loading

a）裸眼井b）套管井

Fig.5Crack propagation process of horizontal wells

3.3單射孔水平井縱剖面壓裂過程模擬

帶有射孔的套管水平井的縱剖面見圖6.追蹤裂紋擴展，持續(xù)施加水壓以觀測水力壓裂裂紋擴展效果.結(jié)構(gòu)尺寸為20.0 m×20.0 m，射孔深度a=1.0 m，網(wǎng)格尺寸為Δx=0.1 m.結(jié)構(gòu)最外層設定為非破壞區(qū)，即s0→∞.加載水壓p從0開始，到t=50 ms為pmax=24 MPa，時間步長Δt=10-5s.

Fig.6Longitudinal section of horizontal well with single perforation and its hydraulic loading curves

單射孔縱剖面水平井水力壓裂過程見圖7.由圖7可知：射孔裂紋在裂紋擴展后發(fā)生裂紋分叉以及多級分叉現(xiàn)象，表明近場動力學方法可以十分方便且有效地模擬裂紋擴展和分叉問題.a）t=0.40 msb）t=0.55 msc）t=0.65 msd）t=70 mse）t=80 msf）t=85 ms圖 7單射孔縱剖面水平井水力壓裂過程

Fig.7Hydraulic fracture process of longitudinal section of horizontal well with single perforation

3.4多射孔水平井縱剖面壓裂過程模擬

水平井在射孔作業(yè)時往往需在套管井井筒內(nèi)壁上射多個孔眼以進行壓裂.該算例考察多射孔水平井的情況.采用40.0 m×20.0 m結(jié)構(gòu)，3射孔同時壓裂.射孔長度a=1.0 m，射孔間距b=5.0 m，見圖8.網(wǎng)格尺寸Δx=0.1 m，時間步長Δt=10-5s.非破壞區(qū)域與水壓加載方式同前例.

圖 8多射孔水平井縱剖面示意

Fig.8Schematic of longitudinal section of horizontal

well with multiple perforations

模擬得到的裂紋擴展過程見圖9.一個顯著的特點是在兩側(cè)裂紋發(fā)生明顯擴展后，中央裂紋才開始緩慢擴展，并且直到最后；相比于兩側(cè)裂紋，中央裂紋的擴展和分叉也很小，表現(xiàn)為兩側(cè)裂紋對中央裂紋具有一定的抑制作用.這表明由于裂紋間的相互影響，多裂紋的擴展過程更為復雜，與單裂紋的擴展過程具有明顯的差異.該算例的模擬結(jié)果對實際工程具有一定的指導意義，表明高密度的初始裂紋分布并不能得到高密度的裂縫網(wǎng)絡.因而，工程上應合理安排射孔排布，以減少成本并達到更好的壓裂效果.

a）t=50 msb）t=70 msc）t=80 msd）t=100 ms圖 9多射孔縱剖面水平井水力壓裂過程模擬結(jié)果

Fig.9Simulation results of hydraulic fracture process in longitudinal section of horizontal well with multiple perforations

4結(jié)束語

利用近場動力學方法便于模擬裂紋萌生、擴展和分叉的顯著優(yōu)點，發(fā)展頁巖水力壓裂過程的數(shù)值模擬方法.模擬結(jié)果表明該方法不僅能模擬頁巖的水力壓裂過程、預測最終的裂縫網(wǎng)絡拓撲，而且能定性揭示該工藝過程的一些力學機理，如：初始射孔能顯著降低起裂水壓，由于裂紋間的相互影響須合理安排射孔間距等.當然，本文工作僅是定性意義上的二維模擬，對于工程亟需的定量的三維模擬仍需深入研究.參考文獻：

[1]OUCHI H， KATIYAR A， YORK J， et al. A fully coupled porous flow and geomechanics model for fluid driven cracks： a peridynamics approach[J]. Computational Mechanics， 2015， 55（3）： 561576. DOI： 10.1007/s0046601511238.

[2]朱海燕， 趙軍， 鄧金根， 等. 水力壓裂過程中微環(huán)隙的產(chǎn)生和擴展的數(shù)值模擬[J]. 計算機輔助工程， 2013， 22（S2）： 443447. DOI： 10.3969/j.issn.10060871.2013.z2.108.

ZHU H Y， ZHAO J， DENG J G， et al. Numerical simulation on microannulus initiation and propagation during hydraulic fracturing[J]. Computer Aided Engineering， 2013， 22（S2）：443447. DOI： 10.3969/j.issn.10060871.2013.z2.108.

[3]董大忠， 鄒才能， 楊樺， 等. 中國頁巖氣勘探開發(fā)進展與發(fā)展前景[J]. 石油學報， 2012， 33（S1）： 107114.

DONG D Z， ZOU C N， YANG H， et al. Progress and prospects of shale gas exploration and development in China[J]. Acta Petrolei Sinica， 2012， 33（S1）： 107114.

[4]BELYTSCHKO T， BLACK T. Elastic crack growth in finite elements with minimal remeshing[J]. International Journal for Numerical Methods in Engineering， 1999， 45（5）： 601620. DOI： 10.1002/（SICI）10970207（19990620）45：5<601：：AIDNME598>3.0.CO；2S.

[5]MOE··S N， DOLBOW J， BELYTSCHKO T. A finite element method for crack growth without remeshing[J]. International Journal for Numerical Methods in Engineering， 1999， 46（1）： 131150. DOI： 10.1002/（SICI）10970207（19990910）46：1<131：：AIDNME726>3.0.CO；2J.

[6]SILLING S A. Reformulation of elasticity theory for discontinuities and longrange forces[J]. Journal of the Mechanics and Physics of Solids， 2000， 48（1）： 175209. DOI： 10.1016/S00225096（99）000290.

[7]HU W， HA Y D， BORARU F. Peridynamic model for dynamic fracture in unidirectional fiberreinforced composites[J]. Computer Methods in Applied Mechanics and Engineering， 2012， 217220： 247261. DOI： 10.1016/j.cma.2012.01.016.

[8]REN H， ZHUANG X， CAI Y， et al. Dualhorizon peridynamics[J]. International Journal for Numerical Methods in Engineering， 2015， 108（12）： 125. DOI： 10.1002/nme.5257.

[9]BESSA M A， FOSTER J T， BELYTSCHKO T， et al. A meshfree unification： reproducing kernel peridynamics[J]. Computational Mechanics， 2014， 53（6）：12511264. DOI： 10.1007/s004660130969x.

[10]SILLING S A， ASKARI E. A meshfree method based on the peridynamic model of solid mechanics[J]. Computers & Structures， 2005， 83（17）： 15261535. DOI： 10.1016/j.compstruc.2004.11.026.

[11]HA Y D， BORARU F. Studies of dynamic crack propagation and crack branching with peridynamics[J]. International Journal of Fracture， 2010， 162（1）： 229244. DOI： 10.1007/s1070401094424.

[12]SONG J H， AREIAS P M， BELYTSCHKO T， et al. A method for dynamic crack and shear band propagation with phantom nodes[J]. International Journal for Numerical Methods in Engineering， 2006， 67（6）： 868893. DOI： 10.1002/nme.1652.（編輯武曉英）第26卷 第1期2017年2月計 算 機 輔 助 工 程Computer Aided EngineeringVol.26 No.1Feb. 2017