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

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

關(guān)鍵詞： 頁巖氣開采； 水力壓裂； 近場(chǎng)動(dòng)力學(xué)； 裂紋擴(kuò)展； 多裂紋； 裂紋追蹤

中圖分類號(hào)： TE357.1文獻(xiàn)標(biāo)志碼： 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

基金項(xiàng)目： 國家自然科學(xué)基金（11232003，11372066）；中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金（DUT15LK07）；

遼寧省教育廳重點(diǎn)實(shí)驗(yàn)室基礎(chǔ)研究項(xiàng)目（LZ2014002）

作者簡介： 吳凡（1992—），男，吉林遼源人，博士研究生，研究方向?yàn)榻鼒?chǎng)動(dòng)力學(xué)，（Email）251386704@qq.com

通信作者： 李書卉（1979—），女，遼寧沈陽人，博士，工程師，研究方向?yàn)椴牧掀茐姆治雠c實(shí)驗(yàn)、并行計(jì)算技術(shù)等，（Email）shuhuili@dlut.edu.cn0引言

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

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

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

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

1近場(chǎng)動(dòng)力學(xué)理論

1.1連續(xù)模型

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

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

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

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

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

2水力壓裂的近場(chǎng)動(dòng)力學(xué)模型

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

Fig.1Schematic of hydraulic fracture process

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

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

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

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

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

3數(shù)值算例

通過幾個(gè)數(shù)值算例考察建立的近場(chǎng)動(dòng)力學(xué)方法模擬水力壓裂過程的效果.為驗(yàn)證近場(chǎng)動(dòng)力學(xué)模型對(duì)裂紋擴(kuò)展模擬的有效性和準(zhǔn)確性，算例1模擬KalthoffWinkler實(shí)驗(yàn)，算例2考察初始裂紋對(duì)所需壓裂水壓的影響情況，算例3考察單射孔水平井縱剖面的水力壓裂裂紋擴(kuò)展情況，算例4進(jìn)行多射孔水平井縱剖面的模擬.

3.1KalthoffWinkler實(shí)驗(yàn)

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

Fig.2KalthoffWinkler experiment setup

Fig.3KalthoffWinkler experiment simulation result obtained by peridynamics method

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

3.2初始裂紋對(duì)壓裂水壓的影響

該算例考慮水平井的橫剖面，分為有初始裂紋的套管井和無初始裂紋的裸眼井.水平井橫剖面示意圖和水壓加載見圖4.結(jié)構(gòu)尺寸為10.00 m×10.00 m，均勻離散下網(wǎng)格尺寸為Δx=0.02 m，節(jié)點(diǎn)數(shù)為65 000個(gè)，中央處水平井直徑為0.60 m，其中，套管井的射孔長度為0.60 m，裸眼井無射孔，對(duì)比2種井各自所需的起裂水壓.水壓加載見圖4b，增量為1 MPa/ms，0.200 ms后恒定（足夠大）.2種井的壓裂效果云圖見圖5.由此可看出：裸眼井的破裂壓力為13.2 MPa，而套管井的射孔在5.8 MPa時(shí)開始發(fā)生裂紋擴(kuò)展，0.250 ms時(shí)發(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.追蹤裂紋擴(kuò)展，持續(xù)施加水壓以觀測(cè)水力壓裂裂紋擴(kuò)展效果.結(jié)構(gòu)尺寸為20.0 m×20.0 m，射孔深度a=1.0 m，網(wǎng)格尺寸為Δx=0.1 m.結(jié)構(gòu)最外層設(shè)定為非破壞區(qū)，即s0→∞.加載水壓p從0開始，到t=50 ms為pmax=24 MPa，時(shí)間步長Δt=10-5s.

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

單射孔縱剖面水平井水力壓裂過程見圖7.由圖7可知：射孔裂紋在裂紋擴(kuò)展后發(fā)生裂紋分叉以及多級(jí)分叉現(xiàn)象，表明近場(chǎng)動(dòng)力學(xué)方法可以十分方便且有效地模擬裂紋擴(kuò)展和分叉問題.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è)時(shí)往往需在套管井井筒內(nèi)壁上射多個(gè)孔眼以進(jìn)行壓裂.該算例考察多射孔水平井的情況.采用40.0 m×20.0 m結(jié)構(gòu)，3射孔同時(shí)壓裂.射孔長度a=1.0 m，射孔間距b=5.0 m，見圖8.網(wǎng)格尺寸Δx=0.1 m，時(shí)間步長Δt=10-5s.非破壞區(qū)域與水壓加載方式同前例.

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

Fig.8Schematic of longitudinal section of horizontal

well with multiple perforations

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

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é)束語

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

[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)生和擴(kuò)展的數(shù)值模擬[J]. 計(jì)算機(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ā)進(jìn)展與發(fā)展前景[J]. 石油學(xué)報(bào)， 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月計(jì) 算 機(jī) 輔 助 工 程Computer Aided EngineeringVol.26 No.1Feb. 2017