蔡申陽(yáng)，尹宏偉*,李長(zhǎng)圣，賈東，汪偉，陳竹新，魏東濤.南京大學(xué)能源科學(xué)研究院，地球科學(xué)與工程學(xué)院，南京00；.中國(guó)石油勘探開(kāi)發(fā)研究院，北京0008；.中國(guó)石油勘探開(kāi)發(fā)研究院西北分院，蘭州700

?

基于離散元數(shù)值模擬的應(yīng)變分析和裂縫預(yù)測(cè)技術(shù)

蔡申陽(yáng)1，尹宏偉1*,李長(zhǎng)圣1，賈東1，汪偉1，陳竹新2，魏東濤3
1.南京大學(xué)能源科學(xué)研究院，地球科學(xué)與工程學(xué)院，南京210023；2.中國(guó)石油勘探開(kāi)發(fā)研究院，北京100083；3.中國(guó)石油勘探開(kāi)發(fā)研究院西北分院，蘭州730022

摘要：應(yīng)變分析與裂縫預(yù)測(cè)技術(shù)在地學(xué)領(lǐng)域具有重要應(yīng)用意義。離散元方法雖然能有效分析含有大量間斷的問(wèn)題，但目前在地學(xué)領(lǐng)域應(yīng)用較少。文中嘗試使用離散元方法表示符合實(shí)際性質(zhì)的巖石，模擬水平擠壓環(huán)境下滑脫構(gòu)造的形成過(guò)程，并對(duì)變形過(guò)程中的應(yīng)變分布變化與裂縫生成規(guī)律進(jìn)行了分析。結(jié)果表明：在擠壓環(huán)境的滑脫構(gòu)造中，裂縫產(chǎn)生的高峰期先于斷層明顯活動(dòng)期，局部區(qū)域內(nèi)聚集的大量裂縫是產(chǎn)生斷層的誘因；已經(jīng)出現(xiàn)明顯活動(dòng)的斷層中產(chǎn)生的裂縫較少。裂縫集中區(qū)域和應(yīng)變集中區(qū)域相互重疊，裂縫越發(fā)育則應(yīng)變?cè)綇?qiáng)烈。受同一個(gè)斷層影響的裂縫首先在斷面上集中出現(xiàn)，隨后產(chǎn)生在斷面周邊區(qū)域；在受斷層影響的范圍內(nèi)，裂縫距離斷面越遠(yuǎn)則形成時(shí)間越晚。該成果還表明離散元方法在應(yīng)變分析與裂縫預(yù)測(cè)研究中具有巨大潛力。

關(guān)鍵詞：離散元方法；應(yīng)變分析；裂縫預(yù)測(cè)

地體應(yīng)變分析及裂縫預(yù)測(cè)技術(shù)在地學(xué)有關(guān)生產(chǎn)研究領(lǐng)域具有重要意義。例如，何雨單和魏春光（2007）認(rèn)為裂縫型油氣藏蘊(yùn)含著巨大的調(diào)整和挖掘潛力，而還原構(gòu)造應(yīng)力場(chǎng)是預(yù)測(cè)裂縫型油藏的重要手段之一；龍鵬宇等（2011）指出頁(yè)巖氣有利勘探目標(biāo)區(qū)應(yīng)首選那些擁有較高滲透能力或具備可改造條件的泥頁(yè)巖裂縫發(fā)育帶；聶海寬等（2009）論證了裂縫的規(guī)模對(duì)頁(yè)巖氣開(kāi)采有雙重影響，適量的裂縫可以提高頁(yè)巖總含氣量，而過(guò)多裂縫則可能導(dǎo)致天然氣大量散失，因此頁(yè)巖氣出產(chǎn)最佳區(qū)域需要深度、溫度、有機(jī)碳含量和裂縫等因素的良好匹配。Gale等（2007）曾分析了美國(guó)Barnett頁(yè)巖的應(yīng)力場(chǎng)和天然裂縫分布，并提出大型開(kāi)放式裂縫的稀疏分布可能不利于水力壓裂開(kāi)采。此外，Lee和Dan（2005）還研究了包括構(gòu)造裂縫在內(nèi)的諸多滑坡誘發(fā)因素，并在此基礎(chǔ)上建立了概率模型，預(yù)測(cè)準(zhǔn)確度達(dá)到了83.47%，表明研究天然裂縫有助于了解滑坡產(chǎn)生的機(jī)制。

目前，通過(guò)物理手段預(yù)測(cè)構(gòu)造裂縫的方式已較為完善。研究者利用光彈模擬物理實(shí)驗(yàn)還原古構(gòu)造應(yīng)力場(chǎng)來(lái)預(yù)測(cè)裂縫的走向與類(lèi)型（Yin and Zhang，2007；單家增等，2005；馬寅生等，2002）；Abdelmalak等（2012）將沙箱模擬實(shí)驗(yàn)與粒子成像測(cè)速（PIV）方法相結(jié)合，建立二維模型，并分析了兩種不同巖漿侵入方式導(dǎo)致的地表破裂變形形態(tài)差異。近年來(lái)，數(shù)值模擬也被廣泛用于構(gòu)造應(yīng)力應(yīng)變分析和裂縫預(yù)測(cè)中，文世鵬和李德同（1996）給出了構(gòu)造裂縫定量預(yù)測(cè)的有限元數(shù)值模擬技術(shù)基本原理和方法，（陳波和田崇魯，1998）則舉出了相應(yīng)的應(yīng)用實(shí)例。有限元方法不能很好地解決含有大量間斷或位移的問(wèn)題，而離散元方法很好地解決了這項(xiàng)不足。離散元方法最先由Cundall 和Strack（1979）完整提出，最初主要用于研究巖塊的節(jié)理破裂以及粒狀物質(zhì)和土壤力學(xué)問(wèn)題的試驗(yàn)?zāi)M中。二十世紀(jì)末，地質(zhì)學(xué)家開(kāi)始將其應(yīng)用于剪切帶問(wèn)題（Iwashita and Oda,2000; Saltzer and Pollard，1992），斷層及斷層相關(guān)褶皺的研究等地質(zhì)構(gòu)造問(wèn)題的研究中（Finch et al.,2003; Hardy and Finchk,2005; Morgan,1997; Saltzer and Pollard,1992; Strayer et al.,2004; Strayer and Suppe，2002;孟令森等，2007;張潔等，2008a，2008b）。離散元方法在一定程度上突破了物理模擬存在的流變學(xué)和比例化問(wèn)題。與其他連續(xù)體數(shù)值模擬方法相比，離散元方法采用顆粒相互作用來(lái)模擬系統(tǒng)的動(dòng)力學(xué)機(jī)制，因此試驗(yàn)者可以對(duì)系統(tǒng)的運(yùn)動(dòng)演化進(jìn)行模擬和觀測(cè)。另外由于它允許顆粒間較大相對(duì)位移，可以更好地模擬高度形變，所以非常適于研究存在大量間斷（如斷層、節(jié)理、破裂）的問(wèn)題。例如，劉泉聲（2011）利用離散元方法分析裂隙巖體等效滲透系數(shù)，得出了應(yīng)力比的增加使水力耦合效果更明顯的結(jié)論；楊艷等（2012）采用BPM（Bonded-Parti?cle Model）模型很好地再現(xiàn)了裂隙巖體水力劈裂的形態(tài)。司馬軍等（2013）模擬了黏性圓形薄層土試樣在粗糙邊界條件下的產(chǎn)生及擴(kuò)展過(guò)程，其裂縫參數(shù)演化規(guī)律與室內(nèi)試驗(yàn)結(jié)果基本一致。孫萍等（2008）通過(guò)離散元分析了裂縫帶錯(cuò)動(dòng)對(duì)地鐵區(qū)間隧道的影響，得出了不同錯(cuò)距工況下隧道襯砌的變形和應(yīng)力。但是，在模擬和預(yù)測(cè)構(gòu)造裂縫領(lǐng)域，離散元的應(yīng)用仍有所欠缺。

本實(shí)驗(yàn)以離散元方法為基礎(chǔ)，采用含固結(jié)力的顆粒材料，模擬擠壓環(huán)境下的滑脫構(gòu)造的形成過(guò)程，并通過(guò)應(yīng)變計(jì)算和裂縫觀測(cè)來(lái)分析裂縫生成規(guī)律，從而探索預(yù)測(cè)構(gòu)造裂縫的新方法

1　離散元的原理

基于離散元方法的計(jì)算模擬一般可以概括為兩部分：第一步，應(yīng)用不同的力-位移法則（即本構(gòu)模型，如固體晶格模型、線性接觸模型等）計(jì)算顆粒間的接觸力；第二步，對(duì)每個(gè)顆粒應(yīng)用牛頓第二定律，更新其位置。反復(fù)執(zhí)行這兩個(gè)步驟，直到計(jì)算結(jié)束。本文的數(shù)值模擬是在開(kāi)源離散元軟件YADE（Yet Another Dynamic Engine）（Kozicki and Donz,2008,2009; Milauer et al.,2014）的框架下進(jìn)行。YADE作為一個(gè)基于面向?qū)ο蟮碾x散元開(kāi)源軟件，易于擴(kuò)展性，設(shè)計(jì)之初就考慮了與各種其他方法（如FEM）耦合，使得許多學(xué)者開(kāi)始采用其做科學(xué)研究（Scholt and Donz，2013;徐文杰等,2012）。較為完美的程序框架和方便的前后處理，目的是讓科學(xué)研究人員集中于真正的科學(xué)問(wèn)題，而不再花費(fèi)太多時(shí)間在程序接口、數(shù)據(jù)的輸入輸出、網(wǎng)格生成及可視化的編程實(shí)現(xiàn)上（Kozicki and Donz,2008）。

1.1本構(gòu)模型

顆粒a、b相互接觸時(shí)，其作用力可以分解為法向力Fn和切向力Fs，法向力Fn可通過(guò)顆粒的法向疊合量得到，切向力通過(guò)增量的形式計(jì)算得到。

顆粒間的法向力與位移的關(guān)系如圖1所示，可分解為壓縮和拉伸兩部分。

圖1　顆粒間的法向力-位移關(guān)系Fig.1 Normal force-displacement relationships between particles

壓縮時(shí)，F(xiàn)n用下式計(jì)算：Fn= kn?D（1）其中，?D為顆粒a和顆粒b的疊合量，kn為法向剛度，其定義如下：

式中，Eeq為等效體積彈性模量，Ra和Rb分別為顆粒a和顆粒b的半徑。

拉伸時(shí)，當(dāng)拉伸距離較小時(shí)，法向力計(jì)算時(shí)的法向剛度取值與壓縮時(shí)一樣。兩顆粒可承受最大的拉力定義為抗拉強(qiáng)度，用下式計(jì)算：

其中，兩顆粒的接觸面積Aint=πr2，這里r取Ra和Rb中的較小值。

當(dāng)兩顆粒間的拉力超過(guò)抗拉強(qiáng)度時(shí)，顆粒間的連接斷裂。能量釋放，需要修正法向剛度，法向力用下式計(jì)算：

當(dāng)?D＞?Drupture時(shí)，顆粒連接發(fā)生拉伸斷裂，顆粒間的接觸力（法向力和切向力）重置為零。

顆粒間的切向力以增量的形式計(jì)算，定義為（Hart et al.,1988）：

ks為切向剛度，?us為顆粒間的切向位移增量。

顆粒間相互作用力應(yīng)符合莫爾庫(kù)倫法則，見(jiàn)圖2。

圖2　顆粒連接的斷裂準(zhǔn)則Fig.2 Breaking rules of a particle bond

顆粒間的最大剪切力為：

φb為內(nèi)摩擦角，c為粘聚力。如果Fs＞Fs,max，顆粒間的連接斷裂，進(jìn)入純滑動(dòng)狀態(tài)。斷裂的連接（拉斷或者切斷）最大剪切力定義為：

φb為連接斷裂后顆粒間的內(nèi)摩擦角。

另外，為了模擬準(zhǔn)靜態(tài)行為，需要給系統(tǒng)加入局部阻尼來(lái)耗散系統(tǒng)的動(dòng)能（Potyondy and Cundall,2004）。

1.2初始連接的生成

一般的離散元方法（Inc,2012; Milauer et al.,2014; Shiu et al.,2008）中，只有兩個(gè)顆粒真實(shí)接觸時(shí)（兩顆粒圓心間的距離小于兩顆粒半徑之和），才生成初始連接。這里重新定義初始連接生成的范圍，給定一個(gè)連接生成的接觸范圍系數(shù)γint，不僅僅真實(shí)接觸的顆粒間會(huì)建立起連接關(guān)系，在給定的范圍內(nèi)的其他鄰居顆粒也會(huì)產(chǎn)生連接關(guān)系，如圖3所示。

初始連接在計(jì)算開(kāi)始前生成，如果兩顆粒圓心間的距離小于平衡距離Deq=γint(Ra+ Rb)，則兩顆粒產(chǎn)生初始連接（Donze et al.,1997）。這里γint取值不能太大，以防產(chǎn)生跨越顆粒的連接，例如當(dāng)Rmax/Rmin= 2時(shí)，γint≤1.5。

為了避免顆粒間產(chǎn)生自鎖力，公式（1）中，相對(duì)位移修正為：

可知，當(dāng)γint= 1時(shí)，顆粒初始連接就是按照一般的離散元方法生成的。

圖3　接觸范圍系數(shù)對(duì)生成連接的影響示例Fig.3 An example of the in fluence of interaction range coefficient on bonds generation

1.3應(yīng)變的計(jì)算

在離散元模型中，作為離散單元的顆粒無(wú)法變形且不可分割，因此顆粒集合體在空間上是離散的。要研究模型的宏觀應(yīng)變，可以采用三角剖分的方法重新構(gòu)造連續(xù)體。YADE采用Delaunay方法（Catalano et al.,2014），以每個(gè)顆粒圓心作為頂點(diǎn)剖分整個(gè)空間，如圖4所示。

通過(guò)計(jì)算剖分單元的形變，可以得到相應(yīng)的應(yīng)變張量，其表示形式如下：

圖4　Delaunay三角剖分法示例Fig.4 An example of Delaunay triangulation

剖分單元的體應(yīng)變?chǔ)舦即為該張量的第一不變量，其值為應(yīng)變張量對(duì)角線上元素的和：

該值反映了剖分單元等比例變形的程度，負(fù)值代表壓縮，正值代表膨脹，且絕對(duì)值越大說(shuō)明變形程度越高。

在得到體應(yīng)變的同時(shí)可以得出球應(yīng)變張量：

從應(yīng)變張量中減去球應(yīng)變張量從而得到偏應(yīng)變張量：

而剖分單元的扭曲程度可以通過(guò)計(jì)算偏應(yīng)變張量的第二不變量（應(yīng)變偏量）得到：

應(yīng)變偏量I2沒(méi)有單位，其值越大，表示剖分單元的扭曲程度越大。

2　構(gòu)造數(shù)值模擬試驗(yàn)

2.1材料參數(shù)

在合適的參數(shù)下，離散元方法可以很好地模擬巖石的變形和斷裂，但是在開(kāi)始試驗(yàn)前需要調(diào)試微觀參數(shù)，使顆粒集合體的宏觀參數(shù)和實(shí)際巖石的力學(xué)性質(zhì)相吻合。通過(guò)大量的模型試驗(yàn)和結(jié)果對(duì)比，筆者得出了一組擠壓實(shí)驗(yàn)結(jié)果比較符合實(shí)際的模型參數(shù)，在該套參數(shù)下的模擬實(shí)驗(yàn)中不會(huì)出現(xiàn)不符合實(shí)際的滑坡或孔穴（Hughes et al.,2014）。本試驗(yàn)所用參數(shù)見(jiàn)表1，顆粒半徑R均勻分布在區(qū)間內(nèi)，所有顆粒的密度為4800 kg/m3，由于顆粒堆積時(shí)存在孔隙，實(shí)際樣品的密度會(huì)小于單個(gè)顆粒的密度。所有顆粒的楊氏模量Eeq為5 GPa，泊松比v為1/3，內(nèi)摩擦角φ為18°，抗張強(qiáng)度t為4.5 MPa,粘聚力c為5 MPa。接觸范圍系數(shù)γint設(shè)置為1.01。

表1　顆粒參數(shù)設(shè)置Table 1 The setting of particles′parameters

Hughes等人曾利用三軸剪切試驗(yàn)導(dǎo)出樣品宏觀參數(shù)，并分析其力學(xué)性質(zhì)是否符合實(shí)際情況（Hughes et al.,2014）。本試驗(yàn)借鑒了該方法，利用三軸剪切試驗(yàn)導(dǎo)出宏觀參數(shù)，并分析其實(shí)際意義，過(guò)程如下：首先將10 000個(gè)顆粒填充于一個(gè)尺寸為1:2:1的剪切試驗(yàn)容器中，如圖5所示。在Yade軟件中，顆粒的絕對(duì)尺寸對(duì)模擬結(jié)果的影響很?。⊿cholt and Donz,2013），因此可以采用內(nèi)壓縮的方式對(duì)樣品施壓。

隨后，讓容器壁以足夠緩慢的速度壓縮樣品，保證其始終處于準(zhǔn)靜止?fàn)顟B(tài)。5個(gè)試驗(yàn)圍壓分別為5 MPa，15 MPa，70 MPa，50 MPa，75 MPa，每個(gè)圍壓下的樣品都是獨(dú)立隨機(jī)生成的。將導(dǎo)出的數(shù)據(jù)繪制成應(yīng)力應(yīng)變曲線圖，曲線在彈性變形階段的斜率即為樣品宏觀楊氏模量，如圖6。

在此基礎(chǔ)上，根據(jù)莫爾-庫(kù)倫破裂準(zhǔn)則可得到樣品的內(nèi)聚力和內(nèi)摩擦角，如圖7。

樣品宏觀力學(xué)性質(zhì)：楊氏模量15.34 GPa,內(nèi)摩擦角28.06°，內(nèi)聚力0.5 MPa。本實(shí)驗(yàn)的參數(shù)與其他離散元構(gòu)造模擬實(shí)驗(yàn)所用參數(shù)基本一致（Finch et al.,2003; Hardy and Finch,2005; Hughes et al.,2014）。在前人的實(shí)驗(yàn)中，用于模擬巖石的樣品的楊氏模量比實(shí)驗(yàn)室結(jié)果大約低一個(gè)數(shù)量級(jí)，這是因?yàn)樽匀唤缰械膸r石體積越大，就可能含有越多的節(jié)理和小斷層，所以一般情況下巖石體積越大則強(qiáng)度越弱（Bieniawski,1984）。

圖5　三軸剪切試驗(yàn)?zāi)Ｐ褪纠龍DFig.5 The example of a trishear model

圖6　三軸壓縮試驗(yàn)下樣品的應(yīng)力應(yīng)變曲線圖（σ1和σ3分別代表軸向應(yīng)力和圍壓，ε1代表軸向應(yīng)變）Fig.6　The stress-strain curve of the sample in trishear tests （The axial stress and confining stress are represented byσ1and σ3respectively.Axial strain is represented by ε1）

圖7　樣品莫爾應(yīng)力圓分析圖（φ代表內(nèi)摩擦角，S0代表內(nèi)聚力）Fig.7 Mohr’s circle analysis on the sample（Theinternalfriction angle and cohesion are represented by φ and S0respectively）

2.2實(shí)驗(yàn)?zāi)Ｐ图斑^(guò)程

實(shí)驗(yàn)?zāi)Ｐ偷倪吔鐬橐婚L(zhǎng)600 m，高200 m，寬4 m的窄盒，顆粒在空間中隨機(jī)生成后自由沉降在盒底，沉積結(jié)束后僅保留距盒底45 m之內(nèi)的顆粒。顆粒組成的初始結(jié)構(gòu)長(zhǎng)600 m，高45 m，寬4 m，如圖8。整個(gè)模型的平均密度為2 691 kg/m3。

實(shí)驗(yàn)中用不同的顏色區(qū)分出6層顆粒，他們的參數(shù)相同且各向同性。所有邊界都是剛性的，且具有不同的摩擦角，其中側(cè)面邊界的摩擦角為18°，底部邊界摩擦角為0.5°。左側(cè)活動(dòng)邊界以1 m/s的速度向右端擠壓模型，直至縮短量達(dá)到30%實(shí)驗(yàn)停止。實(shí)驗(yàn)中每隔40 000步記錄一次變形數(shù)據(jù)。

圖8　均質(zhì)模型的初始結(jié)構(gòu)圖（各層材料參數(shù)均一致；為了方便表示模型形態(tài)和內(nèi)部變形過(guò)程，相鄰層間填充了不同的顏色，以不同灰度體現(xiàn)）Fig.8　Initial configuration of the homogenous model（All layers consist of particles of the same type.The variation of colors,which is based on the greyscale in this figure,are to reflect the model configuration and the internal deformation processes for convenience）

3　實(shí)驗(yàn)結(jié)果分析

3.1變形及應(yīng)變分析

實(shí)驗(yàn)各階段變形情況如圖9。

圖9　模型變形階段圖（間隔200,000步）（為方便表示變形過(guò)程，不同深度的地層著以不同的顏色，但它們的性質(zhì)相同。斷層面的位置和形態(tài)由橙線標(biāo)注，在橙線右上角則對(duì)該斷層加以注釋?zhuān)奖愫笪挠懻摚〧ig.9　Evolution of the deformation with an interval of 200,000 steps (Layers at different depth are colored differently,but their properties are the same.The location and geometry of fault planes are marked out with orange lines.The annotations are given on the upper right of each line for the convenience of later discussion)

橙色線條為斷面所在位置。在樣品縮短量為10.49%時(shí)，緊鄰擠壓端的60 m處一側(cè)出現(xiàn)了兩條傾角約30°的逆斷層f1與f2，其中靠近擠壓端的斷層f1規(guī)模更大，f2規(guī)模較小。在縮短量達(dá)到20.99%時(shí)，f1與f2活動(dòng)強(qiáng)度減弱，在距擠壓端約110 m處發(fā)育第三條主要逆斷層f3。在擠壓試驗(yàn)接近結(jié)束時(shí)，距擠壓端約120 m處出現(xiàn)第四條的逆斷層f4，傾角約40°。該模擬中，斷層最先從靠近擠壓端的一側(cè)產(chǎn)生，且越靠近推覆前端的斷層越新，該現(xiàn)象與孟令森等（2007）得出的結(jié)論相符。

模型變形的應(yīng)變偏量增量如圖10。由應(yīng)變計(jì)算原理可知，處于活動(dòng)中的斷面的應(yīng)變偏量較大，而停止活動(dòng)的斷面上應(yīng)變偏量為零。模型中最新的斷面始終為應(yīng)變偏量最大之處，且最大應(yīng)變所在深度較淺，而早先形成的斷層的活動(dòng)較為輕微。該現(xiàn)象表明在變形過(guò)程中，如果出現(xiàn)了新的斷層則老斷層基本停止活動(dòng)，這與孟令森等（2007）的結(jié)論相符合。

變形過(guò)程中體應(yīng)變?cè)隽咳鐖D11。強(qiáng)烈的體應(yīng)變集中在活動(dòng)斷面的上盤(pán)，且上盤(pán)淺部的應(yīng)變以體膨脹為主。體應(yīng)變強(qiáng)烈的區(qū)域與應(yīng)變偏量較高的區(qū)域相互重疊。

3.2裂縫發(fā)育情況

從變形過(guò)程中破裂數(shù)量統(tǒng)計(jì)圖（圖12）可見(jiàn)4段生成裂縫的高峰期，分別是第0～20 000步，第125 000到160 000步，第235 000到280 000步和第415 000到450 000步。這4個(gè)區(qū)間依次對(duì)應(yīng)在模型變形過(guò)程中先后發(fā)育的4個(gè)主要斷層。

各階段生成的剪裂縫占總裂縫比率見(jiàn)圖13。從圖中可知，在裂縫產(chǎn)生的峰值處，剪裂縫生成比率達(dá)到最大。整個(gè)變形過(guò)程中生成裂縫的總量以張裂居多。

圖10　應(yīng)變偏量增量圖（間隔200,000步）Fig.10 Strain deviator increments with an interval of 200,000 steps

圖11　體應(yīng)變?cè)隽繄D（間隔200,000步）Fig.11 Volumetric strain increments with an interval of 200,000 steps

第一個(gè)裂縫峰期內(nèi)模型的局部變形情況如圖14。第一列圖片為不同破裂模式分布圖，在10 000步時(shí)，新生成裂縫數(shù)達(dá)到最大；靠近擠壓端一側(cè)受到邊界效應(yīng)的影響，因此產(chǎn)生了許多新的裂縫，這解釋了圖12中第一峰期裂縫數(shù)量遠(yuǎn)多于后續(xù)峰期中裂縫數(shù)量的現(xiàn)象。第二列為固結(jié)力分布圖，在該組圖緊靠擠壓端區(qū)域中，大量的固結(jié)鏈接發(fā)生斷裂，其所反映的情況與破裂模式分布圖相符。第三列反映了模型中應(yīng)變偏量的累積變化情況，可見(jiàn)偏應(yīng)變集中的區(qū)域同樣是裂縫集中產(chǎn)生的區(qū)域，二者所處位置具有明顯的對(duì)應(yīng)關(guān)系，符合實(shí)際情況。第四列為模型體應(yīng)變變化情況，體應(yīng)變較強(qiáng)烈的區(qū)域同樣對(duì)應(yīng)了裂縫集中產(chǎn)生的區(qū)域，且應(yīng)變以體膨脹為主。

在第一個(gè)裂縫峰期結(jié)束時(shí)，該處斷層尚未出現(xiàn)明顯位移。

圖15，16，17分別為第二、三、四個(gè)裂縫峰值間距內(nèi)，模型的局部變形情況、破裂分布和應(yīng)變分析圖。它們反映出如下幾個(gè)共同點(diǎn)：（1）裂縫生成的峰期在斷層產(chǎn)生明顯活動(dòng)之前；（2）裂縫和應(yīng)變最先集中在斷面上，屬于同一個(gè)斷層的裂縫距離斷面距離越遠(yuǎn)則生成時(shí)間越晚；（3）裂縫大量發(fā)育的區(qū)域也是應(yīng)變強(qiáng)烈的區(qū)域，新裂縫數(shù)量越多則應(yīng)變?cè)綇?qiáng)烈；（4）同時(shí)在每組圖中，前一段峰期所對(duì)應(yīng)的斷層依舊有明顯的位移，但在其斷面上生成的裂縫數(shù)量遠(yuǎn)遠(yuǎn)小于正在發(fā)育的新斷層。

圖12　裂縫數(shù)量統(tǒng)計(jì)Fig.12 Statistics of the number of fractures

圖13　模擬各階段剪裂縫占總裂縫比率Fig.13 The ratio of shear fractures to total fractures throughout the simulation

圖14　第一裂縫峰期中模型的破裂和應(yīng)變Fig.14　Fractures and strain distribution in the first peak interval

4　討論

本實(shí)驗(yàn)得出的裂縫密度分布規(guī)律與實(shí)際情況相符。前人對(duì)裂縫密度的研究與測(cè)量工作已經(jīng)較為完善。例如，侯貴廷（1994）詳細(xì)介紹了裂縫面密度的測(cè)量方法，即裂縫密度為裂縫總長(zhǎng)度與測(cè)量面的面積之比。利用該方法，Ju等（2014）對(duì)庫(kù)車(chē)凹陷的一處斷層轉(zhuǎn)折褶皺區(qū)域進(jìn)行了裂縫密度測(cè)量，發(fā)現(xiàn)樞紐處的裂縫密度最大，當(dāng)測(cè)量區(qū)域逐漸遠(yuǎn)離斷層時(shí)，測(cè)得的裂縫密度逐漸下降（圖18）。張慶蓮等（2010）對(duì)新疆巴楚地區(qū)走滑斷裂的裂縫密度進(jìn)行了測(cè)量，同樣揭示了越靠近斷層則裂縫密度越大的規(guī)律。Gundmundsson等（2010）采用單位距離內(nèi)的裂縫個(gè)數(shù)來(lái)表示裂縫密度。他們利用該方法在西挪威的一個(gè)斷層周邊統(tǒng)計(jì)了裂縫密度的分布情況，發(fā)現(xiàn)在一定距離內(nèi)離斷層越遠(yuǎn)的區(qū)域裂縫密度越低。Johri等（2014）對(duì)美國(guó)加利福尼亞州近圣安地列斯斷層的SSC儲(chǔ)層的鉆井?dāng)?shù)據(jù)進(jìn)行了研究，同樣發(fā)現(xiàn)了斷層附近區(qū)域的裂縫密度偏高。

圖15　第二裂縫峰期中模型的破裂和應(yīng)變Fig.15　Fractures and strain distribution in the second peak interval

圖16　第三裂縫峰期中模型的破裂和應(yīng)變Fig.16 Fractures and strain distribution in the third peak interval

本實(shí)驗(yàn)中得到的應(yīng)變規(guī)律亦與實(shí)際情況相符合。此前Johri等（2014）在對(duì)SSC儲(chǔ)層的鉆井資料研究中指出斷層核心區(qū)域的應(yīng)變較大；Smart等（2012）對(duì)斷層相關(guān)褶皺進(jìn)行了有限元模擬和應(yīng)力應(yīng)變分析，其結(jié)果同樣表明應(yīng)變主要集中在樞紐區(qū)域。

但是，本實(shí)驗(yàn)也存在不足之處。由于各個(gè)裂縫數(shù)據(jù)都以離散點(diǎn)的形式記錄，而在緊鄰斷層面的空間中的裂縫數(shù)據(jù)點(diǎn)非常多，因此模型中主要斷層面的位置難以精確界定；此外，該裂縫數(shù)據(jù)記錄方式也導(dǎo)致各裂縫的寬度、延伸方向和延伸長(zhǎng)度都無(wú)法得知。在這種情況下，侯貴廷等人所用的裂縫密度測(cè)量公式就不再適用于本實(shí)驗(yàn)了。如何將離散元模型的裂縫密度數(shù)據(jù)與實(shí)際情況做定量對(duì)比，是今后需要解決的問(wèn)題。

圖17　第四裂縫峰期中模型的破裂和應(yīng)變Fig.17 Fractures and strain distribution in the fourth peak interval

圖18　庫(kù)車(chē)凹陷一斷層轉(zhuǎn)折褶皺裂縫密度統(tǒng)計(jì)圖（最靠近斷層區(qū)域的7號(hào)點(diǎn)對(duì)應(yīng)著最高的裂縫密度，而測(cè)量開(kāi)始的1號(hào)點(diǎn)與10號(hào)點(diǎn)上斷層密度最小；Ju et al.,2014）Fig.18 Statistics of fracture densities at different measuring points on a fault-bend fold at Kuqa Depression （Point 7,the one closet to the fault plane,has the highest density of all; While point 1 and 10,the starting and ending measuring points which are far from the fault,bear much lower fracture density）

5　結(jié)論

本實(shí)驗(yàn)利用離散元方法模擬了擠壓環(huán)境下滑脫構(gòu)造的變形過(guò)程，并進(jìn)行了裂縫分析和應(yīng)變分析，得到了如下結(jié)論：

（1）在擠壓環(huán)境的滑脫構(gòu)造中，裂縫產(chǎn)生的高峰期先于斷層明顯活動(dòng)期，局部區(qū)域內(nèi)聚集的大量斷裂是產(chǎn)生斷層的誘因之一。已經(jīng)出現(xiàn)明顯活動(dòng)的斷層中產(chǎn)生的裂縫較少；

（2）裂縫集中區(qū)域和應(yīng)變集中區(qū)域相互重疊，裂縫越發(fā)育則應(yīng)變?cè)綇?qiáng)烈；

（3）受同一斷層影響的裂縫首先在斷面上集中出現(xiàn)，隨后在斷面周邊區(qū)域產(chǎn)生。在斷層影響范圍內(nèi)，裂縫距離斷面越遠(yuǎn)則形成時(shí)間越晚；

本實(shí)驗(yàn)表明，離散元方法很好地模擬了巖體受擠壓產(chǎn)生形變的過(guò)程，并能夠詳細(xì)反映模型內(nèi)部的應(yīng)變和裂縫分布情況。該結(jié)果也同時(shí)揭示了離散元方法在今后裂縫預(yù)測(cè)研究中的巨大潛力。

參考文獻(xiàn)（References）:

陳波，田崇魯.1998.儲(chǔ)層構(gòu)造裂縫數(shù)值模擬技術(shù)的應(yīng)用實(shí)例[J].石油學(xué)報(bào)，19(4): 50-54.

單家增，張占文，陳紹生，等.2005.大民屯凹陷安福屯潛山帶古構(gòu)造應(yīng)力場(chǎng)與裂縫發(fā)育特征的光彈物理模擬實(shí)驗(yàn)研究[J].石油勘探與開(kāi)發(fā)，31(4): 15-18.

何雨丹，魏春光.2007.裂縫型油氣藏勘探評(píng)價(jià)面臨的挑戰(zhàn)及發(fā)展方向[J].地球物理學(xué)進(jìn)展，22(2): 537-543.

侯貴廷.1994.裂縫的分形分析方法[J].應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào)，2(4): 301-306.

劉泉聲，吳月秀，劉濱.2011.應(yīng)力對(duì)裂隙巖體等效滲透系數(shù)影響的離散元分析[J].巖石力學(xué)與工程學(xué)報(bào)，30(1): 176-183.

龍鵬宇，張金川，唐玄，等.2011.泥頁(yè)巖裂縫發(fā)育特征及其對(duì)頁(yè)巖氣勘探和開(kāi)發(fā)的影響[J].天然氣地球科學(xué)，22(3): 525-532.

馬寅生，張興，Others.2002.黃驊坳陷新生代構(gòu)造應(yīng)力場(chǎng)演化的光彈模擬與石油地質(zhì)條件分析[J].地質(zhì)力學(xué)學(xué)報(bào)，8(3): 219-228.

孟令森，尹宏偉，張潔，等.2007.巖石強(qiáng)度和應(yīng)變速率對(duì)水平擠壓變形影響的離散元模擬[J].巖石學(xué)報(bào)，23(11): 2918-2926.

聶海寬，唐玄，邊瑞康.2009.頁(yè)巖氣成藏控制因素及中國(guó)南方頁(yè)巖氣發(fā)育有利區(qū)預(yù)測(cè)[J].石油學(xué)報(bào)，30(4): 484-491.

司馬軍，蔣明鏡，周創(chuàng)兵.2013.黏性土干縮開(kāi)裂過(guò)程離散元數(shù)值模擬[J].巖土工程學(xué)報(bào)，35(zk2): 286-291.

孫萍，彭建兵，范文.2008.地裂縫錯(cuò)動(dòng)對(duì)地鐵區(qū)間隧道影響的三維離散元分析[J].工程地質(zhì)學(xué)報(bào)，5: 710-714.

文世鵬，李德同.1996.儲(chǔ)層構(gòu)造裂縫數(shù)值模擬技術(shù)[J].石油大學(xué)學(xué)報(bào):自然科學(xué)版，20(5): 17-24.

徐文杰，張海洋，于玉貞.2012.基于離散元的土石混合體大型直剪數(shù)值試驗(yàn)研究[J].顆粒材料計(jì)算力學(xué)研究進(jìn)展.

楊艷，常曉林，周偉，等.2012.裂隙巖體水力劈裂的顆粒離散元數(shù)值模擬[J].四川大學(xué)學(xué)報(bào)（自然科學(xué)版），44(5): 78-85.

張潔，尹宏偉，孟令森，等.2008a.主動(dòng)底辟鹽構(gòu)造的二維離散元模擬[J].地球物理學(xué)進(jìn)展，23(6): 1924-1930.

張潔，尹宏偉，徐士進(jìn).2008b.用離散元方法討論巖石強(qiáng)度對(duì)主動(dòng)底辟鹽構(gòu)造斷層分布模式的影響[J].南京大學(xué)學(xué)報(bào):自然科學(xué)版，44(6): 642-652.

張慶蓮，侯貴廷，潘文慶，等.2010.新疆巴楚地區(qū)走滑斷裂對(duì)碳酸鹽巖構(gòu)造裂縫發(fā)育的控制[J].地質(zhì)通報(bào)，29(8): 1160-1167.

Abdelmalak M M,Mourgues R,Galland O,et al.2012.Fracture mode analysis and related surface deformation during dyke intrusion: Results from 2D experimental modelling [J].Earth Planet.Sc.Lett.,359: 93-105.

Bieniawski Z L A T.1984.Rock mechanics design in mining and tunnelling [M] A.A.Balkema,1-272.

Catalano E,Chareyre B and Barth E Lemy E.2014.Pore-scale modeling of fluid-particles interaction and emerging poromechanical effects [J].Int.J.Numer.Anal.Met.,38(1): 51-71.

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

Donze F V,Bouchez J,Magnier S A.1997.Modeling fractures in rock blasting [J].Int.J.Rock Mech.Min.,34(8): 1153-1163.

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

Gale J F,Reed R M and Holder J.2007.Natural fractures in the Barnett Shale and their importance for hydraulic fracture treatments [J].AAPG Bull.,91(4): 603-622.

Gudmundsson A,Simmenes T H,Larsen B,et al.2010.Effects of internal structure and local stresses on fracture propagation,deflection,and arrest in fault zones [J].J.Struct.Geol.,32(11): 1643-1655.

Hardy S and Finch E.2005.Discrete-element modelling of detachment folding [J].Basin Res.,17(4): 507-520.

Hart R,Cundall P A and Lemos J.1988.Formulation of a three-dimensional distinct element model-Part II.Mechanical calculations for motion and interaction of a system composed of many polyhedral blocks [Z].Elsevier: 117-125.

Hughes A N,Benesh N P and Shaw J H.2014.Factors that control the development of fault-bend versus fault-propagation folds: Insights from mechanical models based on the discrete element method (DEM) [J].J.Struct.Geol.,68: 121-141.

Inc I C G.2012.PFC2D/3D (Particle Flow Code In 2/3 Dimensions),Version 4.0 [G].Minneapolis.

Iwashita K and Oda M.2000.Micro-deformation mechanism of shear banding process based on modified distinct element method [J].Powder Technol.,109(1): 192-205.

Johri M,Zoback M D and Hennings P.2014.A scaling law to characterize fault-damage zones at reservoir depths [J].AAPG Bull.,98(10): 2057-2079.

Ju W,Hou G and Zhang B.2014.Insights into the damage zones in fault-bend folds from geomechanical models and field data [J].Tectonophysics,610: 182-194.

Kozicki J and Donz E F V.2009.Yade-open dem: an open-source software using a discrete element method to simulate granular material [J].Eng.Computation,26(7): 786-805.

Kozicki J and Donz E F V.2008.A new open-source software developed for numerical simulations using discrete modeling methods [J].Comput.Method Appl.M.,197(49): 4429-4443.

Lee S and Dan N T.2005.Probabilistic landslide susceptibility mapping in the Lai Chau province of Vietnam: focus on the relationship between tectonic fractures and landslides [J].Environmental Geology,48(6): 778-787.

Morgan J K.1997.Studying Submarine Accretionary Prisms in a'Numerical Sandbox': Simulations using the Distinct Element Method [Z].707.

Potyondy D O and Cundall P A.2004.A bonded-particle model for rock [J].Int.J.Rock Mech.Min.,41(8): 1329-1364.

Saltzer S D and Pollard D D.1992.Distinct element modeling of structures formed in sedimentary overburden by extensional reactivation of basement normal faults [J].Tectonics,11(1): 165-174.

Scholt E S L and Donz E F E D E.2013.A DEM model for soft and hard rocks: Role of grain interlocking on strength [J].J.Mech.Phys.Solids.,61(2): 352-369.

Scholt E S L and Donz E F E D E.2013.A DEM model for soft and hard rocks: Role of grain interlocking on strength [J].J.Mech.Phys.Solids.,61(2): 352-369.

Shiu W,Donze F and Daudeville L.2008.Compaction process in concrete during missile impact: a DEM analysis [J].Comput.Concrete.,5(4): 329-342.

Smilauer V V S,Catalano E,Chareyre B,et al.2014.Yade Documentation [M].The Yade Project,129-133.

Smart K J,Ferrill D A,Morris A P,et al.2012.Geomechanical modeling of stress and strain evolution during contractional fault-related folding [J].Tectonophysics,576: 171-196.

Strayer L M,Erickson S G and Suppe J.2004.Influence of growth strata on the evolution of fault-related folds-distinct-element models [J].AAPG Memoir,82: 413-437

Strayer L M and Suppe J.2002.Out-of-plane motion of a thrust sheet during along-strike propagation of a thrust ramp: a distinct-element approach [J].J Struct.Geol.,24(4): 637-650.

Yin X and Zhang Q.2007.Evolution of the Structural Stress and Its Function During Oil/Gas Accumulation in Bohai Sea [J].Marine Geology and Quaternary Geology,27(5): 45.

Technologyof Strain Analysisand Fracture Prediction Basedon DEMNumerical Simulation

CAI Shenyang1,YIN Hongwei1*,Li Changsheng1,JIA Dong1,WANG Wei1CHEN Zhuxin2,WEI Dongtao3
1.Insititute of Energy Sciences,School of Earth Sciences and Engineering,Nanjing University,Nanjing 210023,China; 2.Research Institute of Petroleum Exploration & Development,Beijing 100083,China; 3.PetroChina Exploration Development Research Institute (Northwest),Lanzhou 730022,China

Abstract:The application of strain analysis and fracture prediction technology is of great importance in geosciences.However,the discrete element method (DEM),which is suitable for solving problems with discontinuous property,has not been widely used in this field.This paper presents an attempt to construct DEM model of a bulk of rock with realistic properties,followed by a simulation of the evolution of a compressional detachment structure.The strain distribution and fracture formation are analyzed afterwards.As to the detachment structure,the result shows that the interval in which the fracture generation reaches its peak precedes the one in which the fault appears the most active.Additionally,fractures and strain distribution share the same concentrating areas,and the number of fractures is proportional to the intensity of the strain.Also,fractures that are related to the same fault first concentrate on its fault plane,then appear in the vicinity of that plane.It indicates that,within the range of fault influence,the further a fracture is from the fault plane,the later it is generated.Above all,the result has revealed a great application potential of DEM in the field of strain analysis and fracture prediction.

Key words:discrete element method; strain analysis; fracture prediction

Corresponding author:YIN Hongwei,Professor; E-mail: hwyin@nju.edu.cn

*通訊作者：尹宏偉，男，1971年生，教授，博士，主要從事含油氣盆地構(gòu)造與沉積特征的研究與教學(xué)工作; E-mail: hwyin@nju.edu.cn

作者簡(jiǎn)介：蔡申陽(yáng)，男，1992年生，實(shí)驗(yàn)助理; E-mail: sycai@outlook.com

基金項(xiàng)目：江蘇省科技支撐計(jì)劃項(xiàng)目（BE2013115）；國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃（“973”：2012CB214703）；國(guó)家自然科學(xué)基金項(xiàng)目（41272227；41572187）聯(lián)合資助

收稿日期：2015-11-30；修回日期：2016-01-18

DOI:10.16108/j.issn1006-7493.2015235

中圖分類(lèi)號(hào)：P618.13

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

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