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

?

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

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

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

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

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

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

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

1　離散元的原理

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

1.1本構(gòu)模型

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

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

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

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

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

拉伸時，當拉伸距離較小時，法向力計算時的法向剛度取值與壓縮時一樣。兩顆?？沙惺茏畲蟮睦Χx為抗拉強度，用下式計算：

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

當兩顆粒間的拉力超過抗拉強度時，顆粒間的連接斷裂。能量釋放，需要修正法向剛度，法向力用下式計算：

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

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

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

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

圖2　顆粒連接的斷裂準則Fig.2 Breaking rules of a particle bond

顆粒間的最大剪切力為：

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

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

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

1.2初始連接的生成

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

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

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

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

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

1.3應(yīng)變的計算

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

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

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

剖分單元的體應(yīng)變εv即為該張量的第一不變量，其值為應(yīng)變張量對角線上元素的和：

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

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

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

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

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

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

2.1材料參數(shù)

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

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

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

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

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

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

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

圖6　三軸壓縮試驗下樣品的應(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實驗?zāi)Ｐ图斑^程

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

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

圖8　均質(zhì)模型的初始結(jié)構(gòu)圖（各層材料參數(shù)均一致；為了方便表示模型形態(tài)和內(nèi)部變形過程，相鄰層間填充了不同的顏色，以不同灰度體現(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　實驗結(jié)果分析

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

實驗各階段變形情況如圖9。

圖9　模型變形階段圖（間隔200,000步）（為方便表示變形過程，不同深度的地層著以不同的顏色，但它們的性質(zhì)相同。斷層面的位置和形態(tài)由橙線標注，在橙線右上角則對該斷層加以注釋，方便后文討論）Fig.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%時，緊鄰擠壓端的60 m處一側(cè)出現(xiàn)了兩條傾角約30°的逆斷層f1與f2，其中靠近擠壓端的斷層f1規(guī)模更大，f2規(guī)模較小。在縮短量達到20.99%時，f1與f2活動強度減弱，在距擠壓端約110 m處發(fā)育第三條主要逆斷層f3。在擠壓試驗接近結(jié)束時，距擠壓端約120 m處出現(xiàn)第四條的逆斷層f4，傾角約40°。該模擬中，斷層最先從靠近擠壓端的一側(cè)產(chǎn)生，且越靠近推覆前端的斷層越新，該現(xiàn)象與孟令森等（2007）得出的結(jié)論相符。

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

變形過程中體應(yīng)變增量如圖11。強烈的體應(yīng)變集中在活動斷面的上盤，且上盤淺部的應(yīng)變以體膨脹為主。體應(yīng)變強烈的區(qū)域與應(yīng)變偏量較高的區(qū)域相互重疊。

3.2裂縫發(fā)育情況

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

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

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

圖11　體應(yīng)變增量圖（間隔200,000步）Fig.11 Volumetric strain increments with an interval of 200,000 steps

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

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

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

圖12　裂縫數(shù)量統(tǒng)計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　討論

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

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

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

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

圖18　庫車凹陷一斷層轉(zhuǎn)折褶皺裂縫密度統(tǒng)計圖（最靠近斷層區(qū)域的7號點對應(yīng)著最高的裂縫密度，而測量開始的1號點與10號點上斷層密度最小；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é)論

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

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

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

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

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

參考文獻（References）:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

張慶蓮，侯貴廷，潘文慶，等.2010.新疆巴楚地區(qū)走滑斷裂對碳酸鹽巖構(gòu)造裂縫發(fā)育的控制[J].地質(zhì)通報，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)造與沉積特征的研究與教學工作; E-mail: hwyin@nju.edu.cn

作者簡介：蔡申陽，男，1992年生，實驗助理; E-mail: sycai@outlook.com

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

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

DOI:10.16108/j.issn1006-7493.2015235

中圖分類號：P618.13

文獻標識碼：A

文章編號：1006-7493（2016）01-0183-11