顧穎凡，盧毅，劉兵，劉春,3,*.南京大學(xué)地球科學(xué)與工程學(xué)院，南京003 .國(guó)土資源部地裂縫地質(zhì)災(zāi)害重點(diǎn)實(shí)驗(yàn)室，江蘇省地質(zhì)調(diào)查研究院，南京0049 3.南京大學(xué)（蘇州）高新技術(shù)研究院，蘇州53

?

基于離散元法的水力壓裂數(shù)值模擬

顧穎凡1，盧毅2，劉兵1，劉春1,3,*
1.南京大學(xué)地球科學(xué)與工程學(xué)院，南京210023 2.國(guó)土資源部地裂縫地質(zhì)災(zāi)害重點(diǎn)實(shí)驗(yàn)室，江蘇省地質(zhì)調(diào)查研究院，南京210049 3.南京大學(xué)（蘇州）高新技術(shù)研究院，蘇州215123

摘要：不同條件下水壓裂隙的發(fā)展特性對(duì)有效開采頁(yè)巖氣具有重要的指導(dǎo)作用。針對(duì)巖體在微觀上為顆粒和孔隙的結(jié)構(gòu)系統(tǒng)，提出離散元水力壓裂數(shù)值模擬方法，離散元能量轉(zhuǎn)化和能量守恒計(jì)算方法，建立了相應(yīng)的三維離散元模型。采用自主研發(fā)的三維離散元模擬軟件MatDEM(3D)，通過(guò)控制模型的豎向應(yīng)變與顆粒直徑，來(lái)模擬地層中的應(yīng)力與壓裂速率的變化。模擬結(jié)果表明：（1）水力壓裂產(chǎn)生裂隙的數(shù)量和方向受巖石的各向異性，壓力狀態(tài)和變化速率所影響。（2）裂隙在壓縮波傳播時(shí)發(fā)展，當(dāng)水壓力高速增加時(shí)，誘發(fā)的裂隙數(shù)量增多，并且有效能量（斷裂熱）百分比也隨之增加，壓裂作用也變得更明顯。（3）當(dāng)豎向應(yīng)變?yōu)榱銜r(shí)，50%的裂隙呈垂直狀態(tài)，當(dāng)豎向應(yīng)變?yōu)?1×10(-4)時(shí)，裂隙趨于沿著最大壓力方向發(fā)展，豎向裂隙的百分率增大。數(shù)值模擬和能量分析為定量地研究巖石水力壓裂過(guò)程提供了一個(gè)新的方法。

關(guān)鍵詞：水力壓裂；離散元；MatDEM；能量

1　引言

頁(yè)巖氣主要成分是甲烷，是一種高效、清潔、低碳的理想替代能源。合理利用頁(yè)巖氣資源可以有效的減輕中國(guó)目前的環(huán)境問(wèn)題，提高能源開發(fā)利用整體效益。水力壓裂技術(shù)是頁(yè)巖氣開發(fā)的核心技術(shù)之一（唐穎等，2011）。早在1947年，美國(guó)Stanolind油氣公司首次將水力壓裂技術(shù)投入使用，隨后水力壓裂技術(shù)便受到廣泛關(guān)注并迅速發(fā)展（張搏等，2015）。由于中國(guó)頁(yè)巖氣儲(chǔ)層的地質(zhì)構(gòu)造較為復(fù)雜，大規(guī)模的工業(yè)化開采亟需有關(guān)技術(shù)的研究與發(fā)展（孫可明等，2014）?？茖W(xué)全面地了解不同的巖石性質(zhì)和應(yīng)力條件下，水力壓裂作用機(jī)制與壓裂紋發(fā)展特性至關(guān)重要。

數(shù)值模擬方法被認(rèn)為是認(rèn)識(shí)裂縫擴(kuò)展規(guī)律經(jīng)濟(jì)而有效的方法，已成為一種研究水力壓裂技術(shù)的有效途徑。目前，一些研究學(xué)者通過(guò)運(yùn)用計(jì)算機(jī)數(shù)值模擬軟件來(lái)模擬壓裂過(guò)程，取得了較好的成果（連志龍等，2009；彪仿俊等，2011；彭成勇等，2014）。連志龍等（2009）運(yùn)用有限元數(shù)值模擬軟件ABAQUS，模擬了地層中各種復(fù)雜因素對(duì)裂縫擴(kuò)展的影響，分析得到了起裂壓力、裂縫擴(kuò)展長(zhǎng)度和最大縫寬與地層中各參數(shù)的關(guān)系。彪仿俊等（2011）建立了一套三維有限元模型，對(duì)水力壓裂中水平裂縫的起裂和擴(kuò)展過(guò)程進(jìn)行了數(shù)值模擬，研究各種參數(shù)對(duì)水平縫擴(kuò)展的影響。彭成勇等（2014）建立了斜井中三維裂縫發(fā)展形態(tài)的有限元模型，綜合考慮了巖層滲透性和孔隙度對(duì)裂縫發(fā)展的影響，研究斜井微環(huán)隙在水力壓裂過(guò)程中的擴(kuò)展特性。

目前，水力壓裂模擬大多采用基于有限元等的連續(xù)介質(zhì)方法，而儲(chǔ)層巖石是由顆粒和孔隙構(gòu)成的不連續(xù)介質(zhì)，基于連續(xù)介質(zhì)方法較難以真實(shí)地反映微觀結(jié)構(gòu)的對(duì)裂隙發(fā)展的作用。離散元法（Cundall and Strack,1979）是解決不連續(xù)物質(zhì)問(wèn)題一種重要方法。在離散元法中，模型由一系列顆粒堆積構(gòu)成，這些顆粒通過(guò)不同的組合來(lái)表示連續(xù)的或者非連續(xù)的物質(zhì)。離散元法能模擬地質(zhì)體的大變形、斷裂和破碎的動(dòng)態(tài)演變過(guò)程（Hazzard et al.,2000）。此種方法已經(jīng)被運(yùn)用在模擬和解釋各種地質(zhì)現(xiàn)象，包括地質(zhì)體的破壞問(wèn)題和地質(zhì)體的不連續(xù)性問(wèn)題上，例如水力壓裂的研究（Shimizu et al.,2011; Eshiet et al.,2013）。

本文建立了緊密堆積的三維離散元模型，通過(guò)轉(zhuǎn)化公式來(lái)確定模型中顆粒力學(xué)參數(shù)，構(gòu)建起具有特定彈性和破壞性質(zhì)的離散元模型，并研究了離散元模型中的能量轉(zhuǎn)化和能量守恒。在此基礎(chǔ)上，模擬了不同豎向應(yīng)變與加壓速率下裂隙發(fā)展特性，嘗試在水力壓裂裂縫發(fā)展機(jī)制和裂縫發(fā)展模式預(yù)測(cè)方面上提供新的思路。本文為定量地研究頁(yè)巖氣水力壓裂過(guò)程提供了一個(gè)新的方法，并有待于進(jìn)一步的實(shí)踐與驗(yàn)證。

2　三維緊密堆積離散元模型

研究采用三維緊密堆積離散元模型（圖1a），其由一系列相同尺寸和力學(xué)性質(zhì)的顆粒堆積構(gòu)成。顆粒間通過(guò)法向彈簧力(Fn)相互作用（Hardy and Finch,2006; Yin et al.,2009）：

其中Kn是法向剛度；Xn是法向相對(duì)位移（圖1b）；Xb是斷裂位移。顆粒最初是相互連接，并受拉力或壓力的作用（式1a）。當(dāng)兩個(gè)顆粒之間的法向相對(duì)位移Xn超過(guò)斷裂位移Xb，連接斷裂，顆粒的拉力消失。此時(shí)，Xn大于等于零時(shí)，兩個(gè)顆粒之間Fn為零（式1c），連接斷裂。但是，當(dāng)這兩個(gè)顆?；貜?fù)到壓縮狀態(tài)，顆粒間斥力仍然存在（式1b）。

圖1c表示兩個(gè)顆粒在切向方向上通過(guò)一個(gè)可斷裂的彈簧連接，并模擬剪切變形和剪切力。剪切力（FS）由下式?jīng)Q定（Cundall and Strack,1979; Hardy et al,2009）：

其中KS是剪切剛度，XS是顆粒間切向相對(duì)位移。在數(shù)值模擬中，剪切力通過(guò)每個(gè)時(shí)間步的計(jì)算累積得到（Cundall and Strack,1979）。對(duì)于一個(gè)完整的連接來(lái)說(shuō)，最大的剪切力（FSmax）由庫(kù)倫摩擦力來(lái)確定：

其中Fs0是顆粒之間的剪切阻力；μp是顆粒之間的摩擦系數(shù)；Fn是法向力（壓縮為負(fù)）。

當(dāng)外力超過(guò)式（3）中的FSmax時(shí)，兩個(gè)顆粒之間的連接會(huì)斷開。剪切力(FS)的大小被限制為小于或等于在下式中的最大剪切值(FSmax')：

當(dāng)顆粒間連接斷裂，且顆粒間剪切力超過(guò)最大剪切力FSmax'，顆粒會(huì)發(fā)生相對(duì)位移，它們之間的滑動(dòng)摩擦力即為FSmax'。

3　模型中能量轉(zhuǎn)換與能量守恒

巖體水力壓裂過(guò)程伴隨著彈性能、動(dòng)能、斷裂能和熱能等能量的產(chǎn)生和轉(zhuǎn)化。定量地研究這些能量的變化，有利于認(rèn)識(shí)水力壓裂機(jī)制和規(guī)律。離散元模型中的機(jī)械能包括動(dòng)能，彈性勢(shì)能和重力勢(shì)能（Mora and Place,1998; Place and Mora,1999）。在數(shù)值模擬中，機(jī)械能通過(guò)三種方式轉(zhuǎn)化為熱能：黏滯力的作用、顆粒間連接斷裂作用與摩擦作用。這三種作用分別產(chǎn)生粘滯熱，斷裂熱和摩擦熱。離散元法中，使用時(shí)間步算法來(lái)進(jìn)行計(jì)算模擬（Cundall and Strack,1979; Potyondy and Cundall,2004），而熱量在每個(gè)時(shí)間步中計(jì)算累積。

3.1粘滯熱

巖體中的應(yīng)力波（地震波）由于摩擦等因素作用而衰減，機(jī)械能逐步轉(zhuǎn)化成熱能（Hazzard et al.,2000）。在離散元模型中，采用粘滯阻尼來(lái)減弱模型中的震動(dòng)波，避免系統(tǒng)動(dòng)能不斷增強(qiáng)（Mora and Place,1993,1994; Place et al.,2002; Finch et al.,2003）。粘滯力（Fv）定義為：

其中η是粘滯系數(shù)，v是顆粒速度。由于模擬中的時(shí)間步非常小，認(rèn)為在一個(gè)時(shí)間步中顆粒的速度為常數(shù)。粘滯熱（Qv）在粘滯力作用下生成，由以下方程來(lái)計(jì)算：

其中dx為在當(dāng)前時(shí)間步中顆粒的位移。

3.2斷裂熱

當(dāng)一個(gè)完整的連接斷裂時(shí)，顆粒之間法向彈力和切向彈力產(chǎn)生的彈性勢(shì)能將會(huì)轉(zhuǎn)為熱能而消散。在連接斷開的過(guò)程中，產(chǎn)生的熱量等于減少的彈性勢(shì)能。如果顆粒之間的法向力為拉力，那么當(dāng)顆粒之間相互的連接斷開的時(shí)候法向與切向的彈性勢(shì)能Ee減少到0。斷裂熱（Qb）為法向和切向彈簧的彈性勢(shì)能的總和：

當(dāng)兩連接顆粒處于相互擠壓狀態(tài)，顆粒間剪切力達(dá)到FSmax時(shí)（式3），顆粒連接斷裂并相對(duì)滑動(dòng)，粒間剪切力減小到FSmax'（式4）。此時(shí)，斷裂熱等于顆粒間切向彈簧熱能的減少量，并可由下式計(jì)算得到：

其中KS是顆粒之間的剪切剛度。

3.3摩擦熱

顆粒相互滑動(dòng)時(shí)產(chǎn)生摩擦熱（Qf），其定義為平均滑動(dòng)摩擦力和滑動(dòng)距離（dS）的乘積：

其中Fs1和Fs2分別為當(dāng)前時(shí)間步開始時(shí)和結(jié)束時(shí)的滑動(dòng)摩擦力（即式4中的FSmax'）；dS為滑動(dòng)距離。由于顆粒之間的法向力在每一個(gè)時(shí)間點(diǎn)都不同，按照式4中得出Fs1和Fs2是不同的，所以在方程中使用平均摩擦力來(lái)計(jì)算。

4　水力壓裂離散元數(shù)值模擬

4.1建模與數(shù)值模擬

基于上文所述原理，我們自主研發(fā)了大型三維離散元模擬軟件MatDEM3D（Liu et al.,2013）。在本研究中，采用此軟件建立三維球體孔隙模型，并模擬巖石中的水力壓裂過(guò)程。如圖2a所示，立方體模型的邊長(zhǎng)（L）為0.06 m，它由一系列直徑為0.001 m的顆粒緊密堆積連接組成。這種緊密堆積的離散元模型具有特定力學(xué)性質(zhì)，并可以通過(guò)轉(zhuǎn)換公式確定顆粒間的5個(gè)力學(xué)參數(shù)（Liu et al.,2013）。附錄A給出三維緊密堆積模型的顆粒間力學(xué)參數(shù)與模型力學(xué)性質(zhì)之間的轉(zhuǎn)換公式。表1給出轉(zhuǎn)換公式中輸入的巖石力學(xué)參數(shù)，計(jì)算得到的顆粒間的力學(xué)參數(shù)，以及實(shí)測(cè)的巖石模型力學(xué)參數(shù)。

圖2?。╝）離散元模型的中心被水顆粒填充；（b）模擬過(guò)程中能量變化曲線；（c）在0.6×10-5s時(shí)的動(dòng)能分布；（d）在0.6×10-5s(εzz為-1×10-4)時(shí)裂隙分布頂視圖；（e）斷裂熱Fig.2 (a) The center of the discrete element model is filled of water particles; (b) Energy conversion curves during simulation; (c) Kinetic energy of model at time 0.6×10-5s; (d) Top view of the fractures 6×10-5s (εzzis -1×10-4); (e) breaking

表1　水力壓裂模型中力學(xué)參數(shù)與顆粒間參數(shù)（顆粒直徑為0.001 m）Table 1　Mechanical propertiesand inter-particle parameters of hydraulic fracturing model (particle diameter is 0.001 m)

為模擬水壓力作用，在模型的中心劃出球形孔區(qū)域，直徑為0.018 m（0.3·L）?？字谐錆M了斷開連接的松散“水顆?！?，即顆粒之間的摩擦力和剪切剛度全為零（μp=0且Ks=0）。由于水沒(méi)有抗拉強(qiáng)度與抗壓強(qiáng)度，斷裂位移和抗剪力都為零（Xb= 0且Fs0=0）。水顆粒的法向剛度可以由水的體積彈性模量（K）來(lái)估計(jì)，體積彈性模量（K）為2.2 GPa。由于水顆粒剪切剛度需為零，當(dāng)水顆粒模型的泊松比（v）設(shè)為0.2時(shí)，由附錄A計(jì)算出的切向剛度為0，符合設(shè)定的要求。據(jù)公式E=3K·(1-2v)，來(lái)計(jì)算水顆粒集合的等價(jià)楊氏模量（E）。將E與v代入附錄A中公式來(lái)計(jì)算法向剛度。轉(zhuǎn)換公式是針對(duì)完整連接的，而數(shù)值模擬中，水顆粒的初始連接全設(shè)為斷裂，實(shí)際泊松比為0.5?；谝陨贤茖?dǎo)，得到水顆粒間的參數(shù)如表1所示。

模擬中，水壓的影響通過(guò)水顆粒之間的相互運(yùn)動(dòng)和相互作用來(lái)實(shí)現(xiàn)，因此需要大量的顆粒。MatDEM3D采用了創(chuàng)新的技術(shù)，可實(shí)現(xiàn)大規(guī)模離散元模擬。此模型由297 578個(gè)顆粒組成，其中4 056個(gè)顆粒在孔隙區(qū)域，512個(gè)顆粒在區(qū)域中心（圖2a中直徑為0.009 m的球形區(qū)域）。原始模型中的鄰近接觸顆粒相互連接。模型的頂端和底部由光滑的平面組成，側(cè)面四個(gè)邊界是無(wú)側(cè)限的。模型被兩個(gè)平面壓縮，直到豎向應(yīng)變（εzz）為-1×10-4。模擬開始時(shí)，核心顆粒直徑增加1%，核心區(qū)域“水壓”增至27 MPa。由于較大的瞬時(shí)水壓作用，模型中產(chǎn)生相應(yīng)的壓縮波（圖2c），顆粒間連接在壓縮波傳播時(shí)斷裂。

圖2b為壓縮波形成和傳播過(guò)程中能量轉(zhuǎn)換曲線。水壓力在0點(diǎn)產(chǎn)生，彈性勢(shì)能（Ee）值跳躍到0.35 J，隨后迅速轉(zhuǎn)化為動(dòng)能（Ek）。時(shí)間點(diǎn)B1，壓縮波到達(dá)孔隙邊界，隨后發(fā)展到孔隙外部（圖2c）。由于水壓力作用，孔隙區(qū)域在時(shí)間點(diǎn)B2擴(kuò)大且開始破裂。在時(shí)間點(diǎn)B3時(shí)，壓縮波到達(dá)模型的邊界，之后斷裂熱保持不變。圖2d所示孔隙區(qū)域破壞特征。圖2e顯示了模型中的斷裂熱。由于模型側(cè)面是無(wú)約束的，當(dāng)壓力波到達(dá)的時(shí)候，裂隙趨于沿著邊界表面發(fā)展。

4.2不同豎向應(yīng)變與壓裂速率

為研究不同的豎向應(yīng)變與壓裂速率對(duì)裂隙數(shù)量與分布的影響，模擬了不同參數(shù)下的破裂過(guò)程。如圖2a為豎向應(yīng)變?yōu)榱銜r(shí)模型H0的頂視圖，其破裂模式與圖2d非常的相似。在低壓裂速率下，經(jīng)過(guò)50個(gè)時(shí)間步，中心顆粒的直徑增加了1%。圖3b-c分別為低壓裂速率豎向應(yīng)變-1×10-4和0時(shí)的破裂結(jié)果。當(dāng)壓裂速率很低時(shí)，所產(chǎn)生的壓縮波比高速率下弱的多。所以，在圖3b-c中裂隙數(shù)量明顯低于在圖2d和3a中的數(shù)量。

圖3　不同條件下的裂隙分布頂視圖，暗色的部分為豎向發(fā)展裂隙。Fig.3 Top views of fractures under different conditions.Note that dark segments are vertical fractures.

表2　在不同豎向應(yīng)變與壓裂速率下的裂隙數(shù)量與能量Table 2 Fracture numbers and energy of simulations with different vertical strains and hydraulic speeds

表2給出了數(shù)值模擬得到的裂隙數(shù)量和能量值。表格中裂隙數(shù)是根據(jù)顆粒間連接的斷開情況判斷的，模型中兩個(gè)顆粒間連接斷開即視作產(chǎn)生一個(gè)裂隙?？梢钥闯?，模型H1、H0產(chǎn)生的裂隙數(shù)量是L1、L0的5.9倍。表2給出有效能百分比，其定義為斷裂熱與輸入彈性勢(shì)能的比值。由表得，H1、H0的平均有效能百分率為2.71%，大約為L(zhǎng)1、L0的三倍。因此，當(dāng)壓裂速率較高時(shí)，裂隙發(fā)育數(shù)量增多，且能量輸入更為有效。表2說(shuō)明當(dāng)豎向應(yīng)變?yōu)榱銜r(shí)，50%的裂隙呈垂直狀態(tài)，由于沒(méi)有豎向壓應(yīng)變時(shí)，顆粒相對(duì)松散，產(chǎn)生裂隙數(shù)量較多，無(wú)側(cè)限時(shí)，水平方向應(yīng)變大，容易產(chǎn)生豎向的裂隙。當(dāng)豎向應(yīng)變?yōu)?1×10-4時(shí)，裂隙趨于平行最大壓應(yīng)力方向發(fā)展，豎向裂隙的百分率較大。

5　結(jié)論

針對(duì)巖體在微觀上為顆粒和孔隙的結(jié)構(gòu)系統(tǒng)的特點(diǎn)，本文采用自主研發(fā)的三維離散元模擬軟件MatDEM3D完成一系列數(shù)值模擬。模擬過(guò)程中，隨著邊界位移的逐漸增加，模型中的總能量增量也逐步的增加，且總是與外力做的功相等。當(dāng)模型破裂，儲(chǔ)存的應(yīng)力能被釋放后轉(zhuǎn)換成動(dòng)能與熱能。但是，機(jī)械能與熱能的總量總為常數(shù)。模擬結(jié)果顯示，模型是一個(gè)孤立的體系，遵循能量守恒規(guī)律。

在本研究中，通過(guò)水顆粒之間相對(duì)運(yùn)動(dòng)和相互作用來(lái)產(chǎn)生水壓力和壓裂作用。模擬結(jié)果顯示，水力壓裂裂隙的數(shù)量和方向受巖石的各向異性，水壓力大小和水壓力增加速率影響等因素影響。當(dāng)壓裂速率很高時(shí)，裂隙數(shù)量增多，并且有效能百分比也隨之增加，壓裂的作用更明顯。本文的初步研究表明，基于MatDEM3D的離散元水力壓裂模型能有效實(shí)現(xiàn)壓裂過(guò)程動(dòng)態(tài)模擬和能量分析，為頁(yè)巖氣水力壓裂數(shù)值模擬研究提供了一個(gè)新的途徑。相關(guān)理論和模型還有待于深入研究，以及進(jìn)一步實(shí)踐的檢驗(yàn)。

附錄A：巖石力學(xué)性質(zhì)與顆粒間力學(xué)參數(shù)的轉(zhuǎn)換公式。顆粒間的法向剛度（Kn），切向剛度（Ks），斷裂位移（Xb），初始剪切抗剪力（Fs0），摩擦系數(shù)（μp）可以通過(guò)楊氏模量（E），泊松比（v），抗拉強(qiáng)度（Tu），抗壓強(qiáng)度（Cu）和內(nèi)摩擦系數(shù)（μi）得到：

參考文獻(xiàn)(References)：

彪仿俊，劉合，張士誠(chéng)，等.2011.水力壓裂水平裂縫影響參數(shù)的數(shù)值模擬研究[J].工程力學(xué)，28(10): 228-235.

連志龍，張勁，王秀喜，等.2009.水力壓裂擴(kuò)展特性的數(shù)值模擬研究[J].巖土力學(xué)，30(1): 169-174.

彭成勇，朱海燕，劉書杰，等.2014.斜井水力壓裂三維裂縫動(dòng)態(tài)擴(kuò)展數(shù)值模擬[J].科技導(dǎo)報(bào)，32(2): 37-43.

唐穎，唐玄，王廣源，等.2011.頁(yè)巖氣開發(fā)水力壓裂技術(shù)綜述[J].地質(zhì)通報(bào)，30(2/3): 393-399.

孫可明，王松，張樹翠.2014.頁(yè)巖氣儲(chǔ)層水力壓裂裂紋擴(kuò)展數(shù)值模擬[J].遼寧工程技術(shù)大學(xué)學(xué)報(bào)（自然科學(xué)版），33(1): 5-10.

張搏，李曉，王宇，等.2015.油氣藏水力壓裂計(jì)算模擬技術(shù)研究現(xiàn)狀與展望[J].工程地質(zhì)學(xué)報(bào)，23(2): 301-310.

Cundall P and Strack O.1979.A discrete element model for granular assemblies [J] Geotechnique,29(1): 47-65.

Eshiet K,Sheng Y,and Ye J.2013.Microscopic modelling of the hydraulic fracturing process [J].Environmental Earth Sciences,68(4)：1169-1186.

Finch E,Hardy S and Gawthorpe R.2003.Discrete element modelling of contractional fault-propagation folding above rigid basement fault blocks [J].Journal of Structural Geology,25(4): 515-528.

Hardy S and Finch E.2006.Discrete element modelling of the influence of cover strength on basement-involved fault-propagation folding [J].Tectonophysics,415(1-4): 225-238.

Hardy S,McClay K and Munoz J.2009.Deformation and fault activity in space and time in high-resolution numerical models of doubly vergent thrust wedges [J].Marine and Petroleum Geology,26(2): 232-248.

Hazzard J,Young R P and Maxwell S C.2000.Micromechanical modeling of cracking and failure in brittle rocks [J].Journal of Geophysical Research-Solid Earth,105(B7): 16683-16697.

Liu C,Pollard D D and Shi B.2013.Analytical solutions and numerical tests of elastic and failure behaviors of close-packed lattice for brittle rocks and crystals [J].Journal of Geophysical Research-Solid Earth,118: 71-82.

Mora P and Place D.1993.A lattice solid model for the non-linear dynamics of earthquakes [J].International Journal of Modern Physics C-Physics and Computers,4(6): 1059-1074.

Mora P and Place D 1994.Simulation of the frictional stick-slip instability [J].Pure and Applied Geophysics,143 (1-3): 61-87.

Mora P and Place D.1998.Numerical simulation of earthquake faults with gouge: toward a comprehensive explanation for the heat flow paradox [J].Journal of Geophysical Research-SolidEarth,103(B0): 21067-21089.

Place D and Mora P.1999.The lattice solid model to simulate the physics of rocks and earthquakes: Incorporation of friction [J].Journal of Computational Physics,150(2): 332-372.

Place D Lombard F,Mora P and Abe S.2002.Simulation of the microphysics of rocks using LSMearth [J].Pure and Applied Geophysics,159 (9): 1911-1932.

Potyondy D and Cundall P.2004.A bonded-particle model for rock [J].International Journal of Rock Mechanics and Mining Sciences,41(8): 1329-1364.

Shimizu H,Murata S and Ishida T.2011.The distinct element analysis for hydraulic fracturing in hard rock considering fluid viscosity and particle size distribution [J].International Journal of Rock Mechanics and Mining Sciences,48(5): 712-727.

Yin H,Zhang J,Meng L,et al.2009.Discrete element modeling of the faulting in the sedimentary cover above an active salt diapir [J].Journal of Structural Geology,31(9): 989-995.

Numerical Simulationof Hydraulic Fracturing Basedon Discrete Element Method

GU Yingfan1,LU Yi2,LIU Bing1,LIU Chun1,3*
1.School of Earth Sciences and Engineering,Nanjing University,Nanjing 210023,China; 2.Key Laboratory of Earth Fissures Geological Disaster,Ministry of Land and Resource,Geological Survey of Jiangsu Province,Nanjing 210049,China; 3.Nanjing University (Suzhou) High-tech Institute,Suzhou 215123,China

Abstract:The studies on the characteristics of hydraulic fracturing under different conditions have an important meaning to the effective shale gas exploration.As rocks are composed of grains and pores at micro scale,the discrete element method is introduced to simulate hydraulic fracturing processes.The rules of energy calculation and energy conversion in the models are proposed and a 3D discrete element model is built.Furthermore,3D discrete element numerical simulation software“MatDEM(3D)”is developed.The variations of strata stress and hydraulic speed are simulated by varying the model vertical strain and the particle diameters,respectively.The simulation results indicate: (1) the number and direction of hydraulic fractures are influenced by the anisotropic properties of rocks,stress state,and hydraulic speed.(2) Fractures are generated when the compressive wave passes.When hydraulic pressure increases at a high rate,the number of induced fractures is much greater.Furthermore,the percentage of effective energybook=195,ebook=198(fracturing energy) is also greater,and as a result,the efficiency of hydraulic fracturing is higher.(3) When the vertical strain is zero,50% of the fractures are vertical.When the vertical strain is -1×10(-4),the fractures tend to develop along the direction of the maximum compressive stress,and the percentages of vertical fractures are greater.The numerical simulation and energy analysis provide a new method for researching the hydraulic fracturing processes of rocks.

Key words:hydraulic fracturing; discrete element method; MatDEM; energy

*通訊作者：劉春，男，副教授，碩士生導(dǎo)師，長(zhǎng)期從事計(jì)算工程地質(zhì)方面研究；E-mail: chunliu@nju.edu.cn

作者簡(jiǎn)介：顧穎凡，男，碩士研究生，主要從事工程地質(zhì)數(shù)值模擬方面研究；E-mail: guyingfan_nju@163.com

基金項(xiàng)目：國(guó)家自然科學(xué)基金青年項(xiàng)目（41302216）；江蘇省自然科學(xué)基金青年項(xiàng)目（BK20130377）；江蘇省科技支撐計(jì)劃（BE2013115）資助

收稿日期：2015-11-20；修回日期：2016-02-15

DOI:10.16108/j.issn1006-7493.2015229

中圖分類號(hào)：TE357.1

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1006-7493（2016）01-0194-06

Correpsponding author: LIU Chun,Associate Professor,E-mail: chunliu@nju.edu.cn