李長(zhǎng)圣, 尹宏偉, 徐雯嶠, 吳珍云, 管樹(shù)巍, 賈 東, 任 榮

基于離散元的擠壓構(gòu)造定量分析與模擬

李長(zhǎng)圣1, 2, 3, 4, 尹宏偉3*, 徐雯嶠3, 吳珍云1, 2, 管樹(shù)巍4, 賈 東3, 任 榮4

(1. 核資源與環(huán)境國(guó)家重點(diǎn)實(shí)驗(yàn)室, 江西 南昌 330013; 2. 東華理工大學(xué)地球科學(xué)學(xué)院, 江西 南昌 330013; 3. 南京大學(xué) 地球科學(xué)與工程學(xué)院, 江蘇 南京 210023; 4. 中國(guó)石油勘探開(kāi)發(fā)研究院, 北京 100083)

隨著離散元理論和計(jì)算機(jī)技術(shù)的發(fā)展, 離散元方法已經(jīng)廣泛應(yīng)用到不同尺度的構(gòu)造模擬中。相較于傳統(tǒng)的沙箱模擬實(shí)驗(yàn), 離散元可以更精確地控制實(shí)驗(yàn)的邊界條件, 定量的分析構(gòu)造變形過(guò)程, 有助于從細(xì)觀尺度深入認(rèn)識(shí)地層力學(xué)性質(zhì)對(duì)構(gòu)造變形過(guò)程的影響。本文系統(tǒng)闡述了基于離散元的構(gòu)造變形定量分析方法, 結(jié)合一個(gè)典型的擠壓構(gòu)造離散元數(shù)值模擬實(shí)驗(yàn), 模擬了水平擠壓環(huán)境下構(gòu)造的形成過(guò)程, 并對(duì)變形過(guò)程中的應(yīng)力、應(yīng)變分布變化與裂縫生成規(guī)律進(jìn)行分析, 取得以下認(rèn)識(shí): (1) 該模擬實(shí)驗(yàn)的構(gòu)造以前展式的逆沖疊瓦斷層為主, 斷層從擠壓端到遠(yuǎn)離擠壓端依次活動(dòng); (2) 裂縫與斷層形成有密切關(guān)系, 局部區(qū)域內(nèi)聚集的大量裂縫是產(chǎn)生斷層的誘因; (3) 斷層形成初期,斷層內(nèi)物質(zhì)運(yùn)動(dòng)距離較小, 產(chǎn)生裂縫增量最多; 斷層活動(dòng)后期, 裂縫增量減少; (4) 體積應(yīng)變可以表征裂縫類型(拉張或壓縮), 變形應(yīng)變可以區(qū)分正向和反向逆沖斷層; (5) 平均應(yīng)力大小與地形起伏呈正相關(guān)關(guān)系, 最大剪切應(yīng)力持續(xù)在將要形成的新斷層處累積, 直至該新斷層形成, 最大剪切應(yīng)力繼而消散, 繼續(xù)往前傳播, 在下一個(gè)將要形成的新斷層處累積。研究結(jié)果表明離散元方法在構(gòu)造變形、應(yīng)力應(yīng)變與裂縫預(yù)測(cè)定量研究中具有巨大潛力。

離散元; 構(gòu)造變形; 擠壓構(gòu)造; 應(yīng)力應(yīng)變; 定量分析

0 引 言

離散元方法(discrete element method, 簡(jiǎn)稱DEM)是Cundall and Strack (1979)提出的用于研究巖土體變形的一種方法, 其基本思想是把自然界的巖土體看成是離散的單元幾何體。離散元方法可以有效地模擬沉積地層中出現(xiàn)的斷層及斷層相關(guān)褶皺等脆性變形(Dean et al., 2013; 李長(zhǎng)圣, 2019), 廣泛應(yīng)用于解決構(gòu)造相關(guān)的地質(zhì)問(wèn)題。例如, 分析鹽刺穿過(guò)程中沉積蓋層的破裂機(jī)制(張潔等, 2008; Yin et al., 2009), 研究地層內(nèi)聚力對(duì)構(gòu)造形態(tài)及應(yīng)力應(yīng)變分布的影響(Morgan, 2015), 揭示水平擠壓環(huán)境下滑脫層對(duì)構(gòu)造變形過(guò)程中的應(yīng)變分布變化與裂縫生成規(guī)律(蔡申陽(yáng)等, 2016), 解析裂陷盆地演化過(guò)程中分層伸展疊加變形的動(dòng)力學(xué)演化過(guò)程(周易等, 2019), 分析區(qū)域應(yīng)力場(chǎng)作用下反轉(zhuǎn)構(gòu)造特性及正斷層反轉(zhuǎn)控震機(jī)制(吳珍云等, 2019), 揭示純走滑拉分盆地發(fā)育過(guò)程中的斷裂擴(kuò)展和連接過(guò)程(劉源和Konietzky, 2019), 探討庫(kù)車(chē)前陸沖斷帶西部鹽構(gòu)造形成的控制因素及其形成機(jī)理(李維波等, 2017; 李江海等, 2020; 徐雯嶠等, 2020)。

Schreurs et al. (2006, 2016)通過(guò)對(duì)比全球多個(gè)物理模擬實(shí)驗(yàn)室的擠壓構(gòu)造實(shí)驗(yàn), 發(fā)現(xiàn)在相同的實(shí)驗(yàn)條件和方法下, 不同的實(shí)驗(yàn)?zāi)M結(jié)果無(wú)法完全一致。而DEM在相同的初始條件(顆粒位置和半徑相同)和邊界條件下, 實(shí)驗(yàn)完全可重復(fù)。其中, DEM顆粒材料屬性通過(guò)細(xì)觀參數(shù)標(biāo)定, 可選擇的材料較多。并且, 所有的變量(如位移、應(yīng)力、應(yīng)變、速度等信息)都可以實(shí)時(shí)監(jiān)測(cè)。在物理模擬中, 變形需要通過(guò)圖像分析(如粒子圖像測(cè)速)(沈禮等, 2012)、激光掃描(微型激光測(cè)高)(Liu et al., 2013)、利用計(jì)算機(jī)X射線斷層掃描技術(shù)進(jìn)行體掃描(李長(zhǎng)圣, 2014; 李長(zhǎng)圣等, 2014; Li et al., 2016)等方法獲取。如果沒(méi)有這些設(shè)備, 則只能通過(guò)對(duì)模型側(cè)面拍照來(lái)觀察模型側(cè)面形變, 或切割模型觀察其內(nèi)部形變。而DEM則可以給出每個(gè)變形階段的所有信息, 用這些信息可以計(jì)算出系統(tǒng)的應(yīng)力應(yīng)變場(chǎng)。同時(shí), 同一初始模型, 模擬結(jié)果一致, 利于單因素變量(如擠壓速率、古隆起、先存斷層、同構(gòu)造沉積與剝蝕等)研究(李長(zhǎng)圣, 2019; 辛文等, 2020; 徐雯嶠等, 2020; Li et al., 2021; Xu et al., 2021; 鄒瑋等, 2021)。

離散元數(shù)值模擬方法在構(gòu)造形態(tài)、斷裂預(yù)測(cè)及應(yīng)力應(yīng)變的定量分析方面具有明顯的優(yōu)勢(shì)。眾多學(xué)者已經(jīng)將離散元引入構(gòu)造變形的模擬中, 但分析的重點(diǎn)多集中于構(gòu)造幾何形態(tài)的定性解析。本文在闡述離散元原理、模擬過(guò)程的基礎(chǔ)上, 結(jié)合一個(gè)典型的擠壓構(gòu)造模擬實(shí)例, 詳細(xì)剖析了基于離散元的擠壓構(gòu)造實(shí)驗(yàn)過(guò)程及定量分析方法。

1 基本原理

DEM基本思想是將顆粒材料內(nèi)部細(xì)觀尺度的單個(gè)離散顆粒視為一個(gè)離散單元, 通過(guò)一系列離散的單元來(lái)模擬顆粒材料的力學(xué)行為。離散元的求解實(shí)際上是迭代計(jì)算顆粒位移和受力, 可以概括為兩部分: 第一步, 已知的顆粒所受合力和合力矩, 由牛頓第二定律更新每個(gè)顆粒的位置; 第二步, 找到相互接觸的顆粒, 應(yīng)用接觸力學(xué)模型(即力?位移法則)。本文采用Hertz-Mindlin接觸力學(xué)模型, 其在力的計(jì)算方面精確且高效, 同時(shí)引入了粘結(jié)(Bond)接觸力學(xué)模型(簡(jiǎn)稱粘結(jié)模型), 以模型脆性巖石的變形行為, 模擬效果較好。

相互粘結(jié)的兩個(gè)顆粒, 其接觸力的計(jì)算可以分為兩種狀態(tài), 即壓縮與拉伸。當(dāng)兩個(gè)顆粒處于壓縮狀態(tài)下, 接觸力計(jì)算采用Hertz-Mindlin接觸力學(xué)模型; 當(dāng)兩個(gè)顆粒處于拉伸狀態(tài)下, 接觸力計(jì)算采用粘結(jié)模型。

一般地, 規(guī)定兩個(gè)顆粒間壓力為正, 拉力為負(fù)。

式中:U. 顆粒間的法向疊合量;K. 顆粒間的法向接觸剛度。

式中:U. 顆粒間的切向疊合量;K. 顆粒間的切向接觸剛度。

(1) 壓縮狀態(tài)

當(dāng)兩個(gè)顆粒(顆粒O和顆粒A, 圖1)處于壓縮狀態(tài)下時(shí),K由下式計(jì)算,

其中,

式中:O和A分別為顆粒O和A的剪切模量,和分別為顆粒O和A的泊松比。

當(dāng)顆粒處于壓縮狀態(tài)下時(shí),K由下式計(jì)算,

(2) 拉伸狀態(tài)

相互粘結(jié)的兩個(gè)顆粒處于拉伸狀態(tài)下時(shí), 法向接觸剛度K和切向接觸剛度K采用以下公式計(jì)算,

圖1 接觸面定義

式中:E. 粘結(jié)的楊氏模量;G. 粘結(jié)的聚合模量;. 兩個(gè)接觸顆粒的接觸面的面積, 為一個(gè)圓, 該圓半徑簡(jiǎn)化為較小顆粒的半徑o, 如圖1所示,o2;o和A分別為顆粒O和A的半徑。

當(dāng)顆粒處于粘結(jié)狀態(tài)時(shí), 粘結(jié)的抗拉力為,

粘結(jié)的抗剪力為,

式中:T.粘結(jié)的抗拉強(qiáng)度;C.粘結(jié)的剪切強(qiáng)度。相互粘結(jié)的顆粒, 自始至終滿足F≤Fmax并且|F|≤Fmax時(shí), 顆粒處于粘結(jié)狀態(tài)。當(dāng)上述一個(gè)條件不滿足時(shí), 顆粒間的粘結(jié)斷開(kāi), 進(jìn)入非粘結(jié)狀態(tài), 之后不再重新產(chǎn)生粘結(jié)。

當(dāng)兩個(gè)顆粒處于非粘結(jié)狀態(tài)(斷裂)時(shí),Fmax= 0.0, 顆粒間不產(chǎn)生拉力;Fmax=μF, 切向力需滿足|F|≤Fmax。

本文的數(shù)值模擬使用李長(zhǎng)圣開(kāi)發(fā)的離散元軟件(李長(zhǎng)圣等, 2017; Li et al., 2018, 2021; 李長(zhǎng)圣, 2019), 選用Hertz-Mindlin接觸力學(xué)模型和粘結(jié)(Bond)模型, 計(jì)算顆粒間的接觸力(Morgan, 2015; 李長(zhǎng)圣, 2019), 采用蛙跳法更新顆粒位移(Itasca Consulting Group, 2008; 李長(zhǎng)圣, 2019)。

2 應(yīng)力應(yīng)變表征

應(yīng)力應(yīng)變作為描述材料受力和變形的基本力學(xué)量, 在連續(xù)介質(zhì)力學(xué)中被廣泛應(yīng)用, 但由于顆粒材料的離散性, 基于連續(xù)介質(zhì)定義的應(yīng)力應(yīng)變無(wú)法直接描述離散顆粒集合體的本構(gòu)特征。本文采用一種在構(gòu)造模擬中應(yīng)用效果較為理想的等效應(yīng)力和等效應(yīng)變定義方法(劉其鵬, 2010)。如沒(méi)有特殊說(shuō)明, 本文提到的應(yīng)力應(yīng)變即為等效應(yīng)力和等效應(yīng)變。

2.1 應(yīng) 變

連續(xù)介質(zhì)力學(xué)中, 應(yīng)變是局部變形區(qū)域的位移梯度(O’Sullivan, 2011)。DEM中, 可以通過(guò)追蹤顆粒每一時(shí)步的位置, 得到顆粒的累積位移, 并分成水平位移和垂直位移兩部分。采用臨近點(diǎn)插值算法(MathWorks, 2015), 將位移投射到規(guī)則的笛卡爾網(wǎng)格中。構(gòu)造研究中, 地層變形一般較大, 需采用有限應(yīng)變理論進(jìn)行計(jì)算分析(Oertel, 1996; Means, 2012)。

集合體隨時(shí)間發(fā)生變化, 可以分解為體積變化導(dǎo)致的形變(體積應(yīng)變)和剪切變形引起的形變(剪切應(yīng)變), 本文采用Morgan (2015)給出的方法計(jì)算應(yīng)變。把每個(gè)顆粒位移分成水平和垂直兩部分, 其二維的位移梯度張量為D,D=?U/?U(指標(biāo)符號(hào),= 1, 2), 其中U為位移矢量網(wǎng)格。有限應(yīng)變張量的不變量用來(lái)表征變形場(chǎng)在時(shí)間和空間上的變化。有限應(yīng)變的第一不變量I是體積應(yīng)變(Malvern, 1969; Jaeger et al., 2009), 由位移梯度張量得出。F=?x/?X,x和X分別為變形后和變形前的顆粒坐標(biāo)。的行列式形式為,

它表征了集合體體積的變化率。

體積應(yīng)變I計(jì)算方法如下:

變形應(yīng)變用有限應(yīng)變張量的第二不變量II′表征, 其中′=–I/2。II′由的分量按下式計(jì)算,

本文采用體積應(yīng)變I(表征體積變化)和有限應(yīng)變張量的第二不變量II′(表征剪切形變程度)描述計(jì)算結(jié)果。

2.2 應(yīng) 力

考慮一個(gè)由圓形顆粒組成的體積為(包括孔隙)的顆粒集合即表征元, 將顆粒集合視為連續(xù)體, 粒材料的應(yīng)力定義為集合內(nèi)的平均應(yīng)力, 假設(shè)個(gè)接觸作用在個(gè)顆粒上, 則(Morgan, 2015; 李長(zhǎng)圣, 2019)

式中:σ. 最大主應(yīng)力;σ. 最小主應(yīng)力; 最大剪切應(yīng)力為max=Δ/2。

式中:r. 顆粒的半徑大小;. 集合體的顆粒數(shù)量。

3 實(shí)驗(yàn)過(guò)程

3.1 參數(shù)標(biāo)定

與有限元、有限差分等連續(xù)介質(zhì)力學(xué)方法不同, 離散元中采用顆粒的細(xì)觀參數(shù)(如顆粒間的剛度系數(shù)、摩擦系數(shù)等)表征顆粒系統(tǒng)的宏觀力學(xué)性質(zhì)。為了模擬自然界地層的構(gòu)造變形行為, 通過(guò)雙軸實(shí)驗(yàn)得到離散元顆粒的細(xì)觀參數(shù)與離散元顆粒體所表現(xiàn)出的宏觀響應(yīng)間的關(guān)系, 即細(xì)觀參數(shù)(如顆粒間的摩擦系數(shù)、剪切模量)與宏觀參數(shù)(如粘聚力、內(nèi)摩擦角等)間的映射關(guān)系。通過(guò)5組雙軸試驗(yàn), 得到如表2所示的細(xì)觀參數(shù)和宏觀參數(shù)對(duì)應(yīng)關(guān)系, 詳細(xì)的測(cè)試過(guò)程和結(jié)果見(jiàn)李長(zhǎng)圣(2019)。雙軸實(shí)驗(yàn)的應(yīng)力應(yīng)變及剪切強(qiáng)度包絡(luò)線見(jiàn)圖2, 顆粒的細(xì)觀參數(shù)見(jiàn)表1, 接觸的粘結(jié)(Bond)參數(shù)見(jiàn)表2。

3.2 擠壓實(shí)驗(yàn)

(1) 在長(zhǎng)60 km, 高15 km的矩形盒子內(nèi), 由92個(gè)直徑160 m的顆粒分別生成左墻和右墻,由749個(gè)直徑80 m顆粒生成底墻。之后, 在該矩形盒子內(nèi), 按照1∶1的比例隨機(jī)生成直徑為120 m和160 m的兩種大小顆粒。新生成顆粒與之前已經(jīng)生成的顆粒進(jìn)行疊合判斷, 如果新生成的顆粒與已生成的顆粒疊合, 則刪除該顆粒, 重新生成一個(gè)顆粒。當(dāng)顆粒生成疊合判斷的嘗試次數(shù)大于2000次(Itasca Consulting Group, 2008; 李長(zhǎng)圣, 2019), 則停止生成顆粒, 共生成25027個(gè)顆粒。

(2) 設(shè)置顆粒間的細(xì)觀參數(shù), 具體值見(jiàn)表1, 所有顆粒的摩擦系數(shù)設(shè)置為0.0。

(3) 顆粒生成之后, 將在重力作用下做自由落體運(yùn)動(dòng), 同時(shí)在阻尼力的作用下, 約5萬(wàn)步后, 沉積穩(wěn)定。之后, 刪除縱坐標(biāo)5 km以上的顆粒(不刪除墻體), 體系上覆地層重量減小, 體系有微小的回彈膨脹, 繼續(xù)計(jì)算1萬(wàn)步, 保證體系穩(wěn)定。最終, 得到長(zhǎng)60 km, 高5 km的初始模型(圖3), 共18606個(gè)顆粒。將左右邊墻和巖層顆粒摩擦系數(shù)設(shè)為0.3, 巖層顆粒細(xì)觀參數(shù)見(jiàn)表1、2。同時(shí), 設(shè)置顆粒顏色, 以便于區(qū)分巖層。

(4) 底墻和右墻固定, 給定左邊墻2 m/s水平向右的恒定速度, 當(dāng)左邊墻向右移動(dòng)20 km時(shí), 結(jié)束計(jì)算。

4 結(jié)果分析

擠壓實(shí)驗(yàn)的構(gòu)造演化過(guò)程見(jiàn)圖4。隨著擠壓的進(jìn)行, 斷層F1、F2、F3和F4依次在模型收縮1 km、4 km、9 km和16 km后開(kāi)始強(qiáng)烈活動(dòng), 產(chǎn)生明顯的斷距, 最終形成典型的疊瓦狀逆沖斷層帶, 該斷層帶整體呈現(xiàn)出由擠壓端向模型內(nèi)部高度逐漸減薄的楔體構(gòu)造形態(tài), 且楔體構(gòu)造外無(wú)明顯的構(gòu)造變形(圖6e)。受擠壓端邊界效應(yīng)的影響, 在靠近擠壓端產(chǎn)生一條反沖斷層, 與物理模擬結(jié)果一致(Sun et al., 2016; 李長(zhǎng)圣, 2019; Cui et al., 2020; Li, 2021)。

圖2 雙軸實(shí)驗(yàn)結(jié)果

表1 顆粒的細(xì)觀參數(shù)

注：按照1∶1的比例隨機(jī)生成直徑為120 m和160 m的兩種大小顆粒。

表2 顆粒間的粘結(jié)參數(shù)

注：壓縮速度取2.0 m/s, 實(shí)驗(yàn)過(guò)程是準(zhǔn)靜態(tài)過(guò)程, 詳細(xì)調(diào)試過(guò)程及結(jié)果見(jiàn)李長(zhǎng)圣(2019)。

圖3 初始模型

圖4 模型收縮量為1 km (a)、4 km (b)、9 km (c)、16 km (d)、20 km (e)時(shí),楔體的構(gòu)造解釋

4.1 楔體構(gòu)造定量分析

楔體的組成要素包括坡頂、坡趾、坡角和坡面(圖4c)。離散元中, 研究對(duì)象為離散顆粒, 可采用基于網(wǎng)格的搜尋地表法識(shí)別出坡面, 采用最遠(yuǎn)距離法搜尋坡趾。與前人研究(Buiter et al., 2006, 2016; Schreurs et al., 2006; Schreurs, 2016; Sun et al., 2016; Wang et al., 2021)相比, 該方法給出了坡趾嚴(yán)格定義, 適合編程實(shí)現(xiàn)自動(dòng)測(cè)量, 排除了人為因素導(dǎo)致的坡角測(cè)量誤差, 詳細(xì)描述見(jiàn)李長(zhǎng)圣(2019)和Li et al. (2021)。

由楔體的高度(圖5a)和寬度(圖5b)變化可知, 楔體高度前期呈線性增加, 后期趨于穩(wěn)定; 楔體寬度則呈現(xiàn)階梯式的增加, 在斷層形成(收縮量1 km、4 km、9 km和16 km)時(shí), 楔體寬度陡增, 物理模擬和數(shù)值模擬均表明楔體寬度存在這種周期性變化的現(xiàn)象(Bigi et al., 2010; Buiter et al., 2016; Schreurs et al., 2016; Sun et al., 2016; Li et al., 2021; Wang et al., 2021), 并且與斷層形成存在密切關(guān)系。與楔體高度和寬度不同, 楔體坡角較快進(jìn)入一個(gè)穩(wěn)定值(圖5c), 總體呈現(xiàn)波動(dòng)狀態(tài), 這種楔體坡角波動(dòng)可以總結(jié)為兩個(gè)狀態(tài): 首先隨著擠壓進(jìn)行楔體逐漸增厚, 坡角顯著增大, 如模型從1 km擠壓到 4 km; 然后, 新的斷層F2形成, 坡角隨即減小, 即模型從4 km收縮到 6 km, 接著楔體繼續(xù)增厚, 坡角顯著增大, 進(jìn)入下一個(gè)演化循環(huán)。新斷層開(kāi)始活動(dòng), 楔體高度持續(xù)增大, 楔體寬度相對(duì)穩(wěn)定, 坡角減小。應(yīng)力累計(jì)、斷裂持續(xù)增加, 形成斷層, 使得坡角發(fā)生波動(dòng), 這種波動(dòng)值大小可能與地層強(qiáng)度有關(guān)。總體上, 坡角趨于一個(gè)穩(wěn)定值, 符合臨界楔理論(Davis et al., 1983)。臨界楔理論已經(jīng)被很多數(shù)值模擬和物理模擬結(jié)果所驗(yàn)證(Buiter, 2012), 成功解釋了如臺(tái)灣造山帶等活動(dòng)邊緣中的許多地質(zhì)現(xiàn)象(Davis et al., 1983; Dahlen, 1984; Dahlen, 1990; Suppe, 2007)。如果沒(méi)有新的物質(zhì)介入, 楔體將沿著滑脫層滑動(dòng)而不產(chǎn)生內(nèi)部形變。而處于次臨界和超臨界狀態(tài)的楔體是不穩(wěn)定的, 其會(huì)調(diào)整內(nèi)部形變來(lái)達(dá)到穩(wěn)定狀態(tài)。次臨界楔體通過(guò)產(chǎn)生新的逆沖斷層或者使老的斷層活化, 而使坡角增加; 超臨界楔體通過(guò)前陸形變, 使坡角減小(孫闖, 2017)。所以, 楔體是以局部滑動(dòng)和產(chǎn)生新的逆沖或使老斷層活化抬升的形式周期性演化的(Morgan, 2015; Sun et al., 2016; 李長(zhǎng)圣, 2019; Li et al., 2021; Wang et al., 2021)。事實(shí)上, 楔體幾何形態(tài)是隨時(shí)間一直變化的, 臨界角只是一個(gè)近似值(Gutscher et al., 1996, 1998)。本次實(shí)驗(yàn)中, 模型的臨界楔角約18° (圖5c), 與物理模擬中石英砂的擠壓實(shí)驗(yàn)相近(李長(zhǎng)圣, 2019; Li et al., 2021)。

圖5 楔體高度(a)、寬度(b)和坡角(c)隨收縮量的變化

4.2 斷裂解譯

擠壓過(guò)程中, 楔體的構(gòu)造以前展式的逆沖疊瓦斷層為主(圖4), 靠近擠壓端的斷層到遠(yuǎn)離擠壓端的斷層依次活動(dòng)。由于顆粒間設(shè)置了粘結(jié)(Bond)力, 相互粘結(jié)的顆粒間具有粘結(jié)力鏈, 使得楔體具有抗剪強(qiáng)度和抗拉強(qiáng)度, 斷層形成時(shí), 粘結(jié)力鏈會(huì)發(fā)生剪切或者拉伸斷裂。

收縮量為1 km、4 km、9 km、16 km、20 km時(shí), 粘結(jié)力鏈分布圖見(jiàn)圖6。模型初始化階段, 模型中相互接觸的顆粒生成粘結(jié)力鏈, 粘結(jié)力鏈數(shù)最大。隨著模型收縮, 粘結(jié)力鏈達(dá)到設(shè)置的抗剪或者抗拉強(qiáng)度, 將有粘結(jié)力鏈斷開(kāi)。顆粒斷開(kāi)后, 可能再次接觸, 但是它們之間不再產(chǎn)生粘結(jié), 這時(shí)它們?yōu)闊o(wú)粘結(jié)接觸。隨著擠壓進(jìn)行, 粘結(jié)力鏈逐漸斷開(kāi), 依次形成斷層F1、F2、F3和F4, 這些斷層處, 地層裂縫發(fā)育密集, 斷層間地層裂縫相對(duì)稀疏。模型收縮量達(dá)到20 km時(shí), F4斷層前端形成一對(duì)“V”字型正反逆沖斷層, 表明斷層并不是從底墻某一破裂點(diǎn)發(fā)育而來(lái), 而是從淺地表發(fā)育密集小斷裂, 并逐漸先下延伸貫穿為斷層, 與Wang et al. (2021)研究結(jié)果不同。

圖6 收縮量為1 km (a)、4 km (b)、9 km (c)、16 km (d)、20 km (e)時(shí), 粘結(jié)力鏈分布(藍(lán)色為粘結(jié)力鏈, 黃色為無(wú)粘結(jié)的接觸)

每擠壓1 km為一個(gè)粘結(jié)力鏈數(shù)統(tǒng)計(jì)單位, 每擠壓1 km, 新斷開(kāi)的粘結(jié)力鏈數(shù)見(jiàn)圖7。例如, 模型從0 km收縮到1 km, 新斷開(kāi)的粘結(jié)力鏈數(shù)為1611個(gè), 從3 km到4 km新斷開(kāi)的粘結(jié)力鏈數(shù)為1826個(gè)。模型從0 km收縮到1 km, 斷層F1開(kāi)始活動(dòng); 從3 km收縮到4 km, 斷層F2開(kāi)始活動(dòng); 從8 km收縮到9 km, 斷層F3開(kāi)始活動(dòng); 從15 km收縮到16 km, 斷層F4開(kāi)始活動(dòng)。新斷層的形成具有瞬時(shí)性, 以圖7中F3斷層為例, 新斷層(F3)開(kāi)始活動(dòng)時(shí), 新斷開(kāi)的粘結(jié)力鏈最大, 之后, 新斷開(kāi)的粘結(jié)力鏈數(shù)逐漸減少。也就是說(shuō), 隨著模型逐漸收縮, F3斷層處應(yīng)力累積導(dǎo)致F3形成, 形成密集地層裂縫, 釋放大量能量。之后, 該斷層將逐漸緩慢釋放能量, 直到更新的斷層F4形成, 此時(shí)老斷層F3停止活動(dòng)。

模型收縮過(guò)程中, 斷層的滑移曲線見(jiàn)圖8。隨著模型收縮, 斷層F1、F2、F3和F4依次產(chǎn)生, 新斷層產(chǎn)生, 老斷層活動(dòng)逐漸停止。模型在收縮到3 km時(shí), 斷層F2可能已經(jīng)開(kāi)始活動(dòng), 模型收縮到4 km時(shí), F2斷層在變形、力鏈和應(yīng)變圖上可見(jiàn); F2斷層形成時(shí), F1斷層仍在活動(dòng), 相鄰斷層活動(dòng)時(shí)間有疊合(圖8灰色部分), 新斷層產(chǎn)生, 老斷層逐漸趨于穩(wěn)定。F3和F2、F4和F3活動(dòng)時(shí)間也有同樣的現(xiàn)象。

選取斷層F1、F2、F3和F4處的四個(gè)監(jiān)測(cè)點(diǎn)PF1、PF2、PF3和PF4, 繪制了它們的運(yùn)動(dòng)路徑(圖9a)。隨著模型收縮, 點(diǎn)PF1累計(jì)位移量最大, PF2、PF3和PF4依次減小(圖9b)。為了進(jìn)一步分析斷層內(nèi)粘結(jié)力鏈演化過(guò)程, 統(tǒng)計(jì)了各個(gè)監(jiān)測(cè)點(diǎn)周?chē)? km范圍內(nèi)的粘結(jié)力鏈斷開(kāi)情況(圖9c)。各個(gè)斷層活動(dòng)時(shí), 其內(nèi)部地層發(fā)育大量斷裂縫, 新斷層活動(dòng), 老斷層內(nèi)部斷裂縫分布基本不再變化, 斷層間具有較強(qiáng)的相互獨(dú)立性, 相互間的影響甚微。處于活動(dòng)中的斷層內(nèi), 新斷開(kāi)的粘結(jié)力鏈數(shù)、斷裂數(shù)多, 而停止活動(dòng)的斷層內(nèi)新斷開(kāi)的粘結(jié)力鏈則趨于零。該現(xiàn)象表明在變形過(guò)程中, 如果出現(xiàn)了新的斷層則老斷層基本停止活動(dòng)(孟令森等, 2007; 蔡申陽(yáng)等, 2016)。

4.3 應(yīng)力和應(yīng)變

收縮量不同時(shí)的累積體積應(yīng)變見(jiàn)圖10。由于左墻擠壓作用, 模型整體呈現(xiàn)體壓縮狀態(tài), 靠近擠壓端應(yīng)變較大, 遠(yuǎn)離擠壓端應(yīng)變較小, 斷層內(nèi)部體膨脹和體壓縮是共存的, 楔體淺部呈現(xiàn)體膨脹狀態(tài), 與前人研究結(jié)果一致(Morgan, 2015)。由于深部圍壓大, 模型體積膨脹受限, 深部體膨脹較少。

圖7 每擠壓1 km過(guò)程中, 新斷開(kāi)的粘結(jié)力鏈數(shù)

收縮量不同時(shí), 模型內(nèi)的累積變形應(yīng)變見(jiàn)圖11。變形應(yīng)變強(qiáng)的區(qū)域與體應(yīng)變較高的區(qū)域相互重疊(蔡申陽(yáng)等, 2016)。正向逆沖(紅色)往往伴隨有反向逆沖(藍(lán)色)斷層, 尤其當(dāng)?shù)撞炕搶訌?qiáng)度較弱時(shí), 這種共軛的“X”狀發(fā)育的斷層更為明顯(Morgan, 2015;李長(zhǎng)圣和尹宏偉, 2017; 李長(zhǎng)圣, 2019)。斷層F3(圖11c)和斷層F4(圖11d)的淺部變形應(yīng)變較大, 也表明斷層是從淺部開(kāi)始發(fā)育, 并逐步向深部發(fā)展。

圖8 斷層的滑移曲線隨收縮量的變化(斷層F1、 F2、 F3和F4見(jiàn)圖4, 斷層滑移值選取綠色層測(cè)量)

平均應(yīng)力演化過(guò)程見(jiàn)圖12, 平均應(yīng)力最大的地方均集中在靠近擠壓端的深部地層, 越遠(yuǎn)離擠壓端, 平均應(yīng)力越小, 地層越厚的地方平均應(yīng)力越大。收縮量為20 km時(shí), V1、V2、V3和V4處的平均應(yīng)力隨深度變化曲線見(jiàn)圖13。平均應(yīng)力與自重應(yīng)力分布基本一致, 地層厚度大的地方, 平均應(yīng)力大。選取水平方向H1測(cè)線, 可見(jiàn)平均應(yīng)力分布與地形起伏呈正相關(guān)關(guān)系(圖14)。

隨著地層持續(xù)收縮, 最大剪切應(yīng)力首先在靠近擠壓端深部集中, 形成斷層F1, 之后不斷往擠壓前端傳播的, 持續(xù)在將要形成的新斷層F3處累積(圖15b), 直至斷層F3形成(圖15c), 剪切應(yīng)力開(kāi)始消散(圖15d), 繼續(xù)往前傳播, 在下一個(gè)將要形成的新斷層處累積。最大剪切應(yīng)力是一個(gè)體系的瞬時(shí)狀態(tài), 斷層形成后最大剪切應(yīng)力釋放, 斷層附近最大剪切應(yīng)力集中現(xiàn)象將消失。隨著擠壓進(jìn)行, 最大剪切應(yīng)力隨楔體逐漸向前傳播, 在推覆前端累積, 斷層形成, 應(yīng)力消散, 并在推覆前端繼續(xù)累積。盡管深部最大剪切應(yīng)力較大, 但是變形往往從淺部開(kāi)始(圖11),并逐步向深部擴(kuò)展。

圖a中4個(gè)點(diǎn)分別選取自綠色標(biāo)志層中的斷層F1、F2、F3和F4內(nèi)部的4個(gè)顆粒, 為便于識(shí)別, 這4個(gè)顆粒半徑放大3倍, 藍(lán)色力鏈為監(jiān)測(cè)點(diǎn)周?chē)? km范圍內(nèi)的粘結(jié)力鏈。

藍(lán)色表示體壓縮, 紅色表示體膨脹, 顏色深淺表征體積變化的大小。其中, (e)中半透明黑色為標(biāo)志層。

藍(lán)色表示逆時(shí)針剪切, 代表反沖斷層帶; 紅色表示順時(shí)針剪切, 代表逆沖斷層帶。顏色深淺表征了剪切應(yīng)變的大小。其中, (e)中黑色為標(biāo)志層。

黑色線條為應(yīng)力等值線, 黑色點(diǎn)為|變形應(yīng)變|>4的點(diǎn), 表征了斷層的位置。子圖(e)中, V1、V2、V3和V4為應(yīng)力監(jiān)測(cè)線, 黃色點(diǎn)為PF1、PF2、PF3和PF4位置, 綠色層為標(biāo)志層。

自重應(yīng)力按照ρgh計(jì)算, 其中, 巖層密度ρ取2500 km/m3, 重力加速度g取10, h為地層深度。監(jiān)測(cè)線V1、V2、V3、V4的位置見(jiàn)圖12, 黃色點(diǎn)為PF1、PF2、PF3和PF4的應(yīng)力值。

5 結(jié) 論

近年來(lái), 離散元方法越來(lái)越多地應(yīng)用到構(gòu)造模擬中。本文簡(jiǎn)述了離散元的基本原理, 通過(guò)雙軸實(shí)驗(yàn)標(biāo)定了顆粒的細(xì)觀參數(shù)對(duì)應(yīng)的地層宏觀力學(xué)參數(shù)(粘聚力10.5 MPa, 內(nèi)摩擦角18.6°), 實(shí)現(xiàn)了一個(gè)典型的擠壓構(gòu)造離散元數(shù)值模擬實(shí)驗(yàn), 并給出了基于離散元的構(gòu)造變形定量分析方法, 取得以下認(rèn)識(shí):

(1) 該模擬實(shí)驗(yàn)的構(gòu)造以前展式的逆沖疊瓦斷層為主, 靠近擠壓端的斷層到遠(yuǎn)離擠壓端的斷層依次活動(dòng);

(2) 裂縫與斷層形成有密切關(guān)系, 局部區(qū)域內(nèi)聚集的大量裂縫是產(chǎn)生斷層的誘因;

(3) 斷層形成是瞬時(shí)的, 斷層形成初期, 粘結(jié)力鏈增量最大, 呈現(xiàn)開(kāi)口向上的弧形, 之后, 該新斷層將持續(xù)活動(dòng), 老斷層活動(dòng)基本停止, 新斷裂的粘結(jié)力鏈數(shù)逐漸減小; 斷層之間相互獨(dú)立, 斷層形成初期, 斷層內(nèi)物質(zhì)運(yùn)動(dòng)距離較小, 但是新產(chǎn)生裂縫數(shù)量最多; 斷層活動(dòng)后期, 斷層內(nèi)物質(zhì)隨著模型收縮, 位移呈現(xiàn)線形增加, 新產(chǎn)生的裂縫數(shù)較少。微斷裂往往從淺部發(fā)育并逐步向深部擴(kuò)展, 最終貫穿形成斷層;

圖14 收縮量為20 km時(shí),水平方向平均應(yīng)力分布

黑色線條為應(yīng)力等值線, 黑色點(diǎn)為|變形應(yīng)變|>4的點(diǎn), 表征了斷層的位置。圖(e)中, 綠色層為標(biāo)志層。

(4) 體積應(yīng)變可以表征裂縫類型(拉張或壓縮), 變形應(yīng)變可以區(qū)分正向和反向逆沖斷層;

(5) 隨著地層收縮, 平均應(yīng)力大小與地形起伏呈正相關(guān)關(guān)系, 呈現(xiàn)楔狀體分布, 靠近擠壓端地層厚度最大; 隨著地層持續(xù)收縮, 最大剪切應(yīng)力在遠(yuǎn)離擠壓端逐漸累積, 最大剪切應(yīng)力不斷往擠壓前端傳播的, 持續(xù)在將要形成的新斷層處累積, 直至該新斷層形成, 剪切應(yīng)力開(kāi)始消散, 繼續(xù)往前傳播, 在下一個(gè)將要形成的新斷層處累積。

離散元方法將顆粒集合體模型視為若干離散單元的集合, 允許顆粒間產(chǎn)生大位移, 特別適合用于構(gòu)造變形定量研究。離散元方法在構(gòu)造變形、應(yīng)力應(yīng)變與裂縫預(yù)測(cè)定量研究中具有巨大潛力, 是未來(lái)的構(gòu)造變形定量研究的主要方向之一。

致謝: 本文的數(shù)值計(jì)算在南京大學(xué)高性能計(jì)算中心的計(jì)算集群上完成, 數(shù)值模擬實(shí)驗(yàn)使用課題組自主研發(fā)的離散元數(shù)值模擬軟件完成, 更多構(gòu)造模擬示例見(jiàn)https: //geovbox.com。文中采用的應(yīng)變計(jì)算代碼, 修改自萊斯大學(xué)Julia K. Morgan和Thomas Fournier的腳本, 在此表示感謝。楔體要素測(cè)量程序源碼見(jiàn)https: //github.com/demsheng/Quantitative-wedge。中國(guó)石油大學(xué)(北京)余一欣副教授和劉志娜副教授對(duì)本文進(jìn)行了認(rèn)真而專業(yè)的審閱, 提出了建設(shè)性修改建議, 在此一并致以特別感謝。

蔡申陽(yáng), 尹宏偉, 李長(zhǎng)圣, 賈東, 汪偉, 陳竹新, 魏東濤. 2016. 基于離散元數(shù)值模擬的應(yīng)變分析和裂縫預(yù)測(cè)技術(shù). 高校地質(zhì)學(xué)報(bào), 22(1): 183–193.

李江海, 章雨, 王洪浩, 王殿舉. 2020. 庫(kù)車(chē)前陸沖斷帶西部古近系鹽構(gòu)造三維離散元數(shù)值模擬. 石油勘探與開(kāi)發(fā), 47(1): 65–76.

李維波, 李江海, 王洪浩, 黃少英, 能源. 2017. 庫(kù)車(chē)前陸沖斷帶克拉蘇構(gòu)造帶變形影響因素分析——基于離散元數(shù)值模擬研究. 大地構(gòu)造與成礦學(xué), 41(6): 1001– 1010.

李長(zhǎng)圣. 2014. 含礫滑帶土三軸剪切的精細(xì)數(shù)值模擬研究. 南京: 南京大學(xué)碩士學(xué)位論文.

李長(zhǎng)圣. 2019. 基于離散元的褶皺沖斷帶構(gòu)造變形定量分析與模擬. 南京: 南京大學(xué)博士學(xué)位論文.

李長(zhǎng)圣, 尹宏偉. 2017. 滑脫層強(qiáng)度對(duì)擠壓構(gòu)造的影響: 離散元數(shù)值模擬 // 2017中國(guó)地球科學(xué)聯(lián)合學(xué)術(shù)年會(huì)論文集. 北京: 中國(guó)和平音像電子出版社: 3879–3882.

李長(zhǎng)圣, 尹宏偉, 劉春, 蔡申陽(yáng). 2017. 共享內(nèi)存式并行離散元程序的設(shè)計(jì)與測(cè)試. 南京大學(xué)學(xué)報(bào)(自然科學(xué)版), 53(6): 1161–1170.

李長(zhǎng)圣, 張丹, 王宏憲, 獨(dú)莎莎. 2014. 基于CT掃描的土石混合體三維數(shù)值網(wǎng)格的建立. 巖土力學(xué), 35(9): 2731–2736.

劉其鵬. 2010. 基于平均場(chǎng)理論的顆粒材料離散顆粒集合- Cosserat連續(xù)體模型多尺度模擬. 大連: 大連理工大學(xué)博士學(xué)位論文.

劉源, Konietzky H. 2019. 基于離散元方法對(duì)走滑拉分盆地演化及次級(jí)斷裂擴(kuò)展過(guò)程研究. 地質(zhì)力學(xué)學(xué)報(bào), 25(5): 840–852.

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

沈禮, 賈東, 尹宏偉, 孫闖, 張勇, 范小根. 2012. 基于粒子成像測(cè)速(PIV)技術(shù)的褶皺沖斷構(gòu)造物理模擬. 地質(zhì)論評(píng), 58(3): 471–480.

孫闖. 2017. 龍門(mén)山褶皺沖斷帶構(gòu)造物理模擬研究. 南京: 南京大學(xué)博士學(xué)位論文.

吳珍云, 尹宏偉, 李長(zhǎng)圣, 孫業(yè)君, 李麗梅, 杜航, 何奕成, 劉博雅. 2019. 斷陷盆地正反轉(zhuǎn)構(gòu)造實(shí)驗(yàn)?zāi)M及對(duì)茅山斷裂帶的啟示. 南京大學(xué)學(xué)報(bào)(自然科學(xué)版), 55(5): 869–878.

辛文, 陳漢林, 安凱旋, 張欲清, 楊樹(shù)鋒, 程曉敢, 林秀斌. 2020. 基于離散元數(shù)值模擬的西南天山山前沖斷帶構(gòu)造變形控制因素研究. 地質(zhì)學(xué)報(bào), 94(6): 1704–1715.

徐雯嶠, 汪偉, 尹宏偉, 賈東, 李長(zhǎng)圣, 楊庚兄, 李剛. 2020. 庫(kù)車(chē)坳陷東西段鹽下構(gòu)造變形差異演化數(shù)值模擬分析. 地質(zhì)學(xué)報(bào), 94(6): 1740–1751.

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

周易, 于福生, 劉志娜. 2019. 分層伸展疊加變形二維離散元模擬. 大地構(gòu)造與成礦學(xué), 43(2): 213–225.

鄒瑋, 余一欣, 劉金水, 蔣一鳴, 唐賢君, 陳石, 余浪. 2021. 東海盆地西湖凹陷中央反轉(zhuǎn)構(gòu)造帶發(fā)育主控因素及寧波背斜形成過(guò)程. 石油學(xué)報(bào), 42(2): 176–185.

Bigi S, Di Paolo L, Vadacca L, Gambardella G. 2010. Load and unload as interference factors on cyclical behavior and kinematics of Coulomb wedges: Insights from sandboxexperiments., 32(1): 28–44.

Buiter S J. 2012. A review of brittle compressional wedge models., 530: 1–17.

Buiter S J, Babeyko A Y, Ellis S, Gerya T V, Kaus B J, Kellner A, Schreurs G, Yamada Y. 2006. The numerical sandbox: Comparison of model results for a shortening and an extension experiment.,,, 253: 29–64.

Buiter S J, Schreurs G, Albertz M, Gerya T V, Kaus B, Landry W, Le Pourhiet L, Mishin Y, Egholm D L, Cooke M. 2016. Benchmarking numerical models of brittle thrust wedges., 92: 140–177.

Cui J, Jia D, Yin H W, Chen Z X, Li Y Q, Wang M M, Fan X G, Shen L, Sun C, Li Z G, Ma D L, Zhang Y K. 2020. The influence of a weak upper ductile detachment on the Longmen Shan fold-and-thrust belt (eastern margin of the Tibetan Plateau): Insights from sandbox experiments., 198, 104220.

Cundall P A, Strack O D. 1979. A discrete numerical model for granular assemblies., 29(1): 47–65.

Dahlen F A. 1984. Noncohesive critical Coulomb wedges: An exact solution.:, 89(B12): 10125–10133.

Dahlen F A. 1990. Critical taper model of fold-and-thrust belts and accretionary wedges., 18(1): 55–99.

Davis D, Suppe J, Dahlen F A. 1983. Mechanics of fold-and- thrust belts and accretionary wedges.:, 88(B2): 1153–1172.

Dean S L, Morgan J K, Fournier T. 2013. Geometries of frontal fold and thrust belts: Insights from discrete element simulations., 53: 43–53.

Gutscher M, Kukowski N, Malavieille J, Lallemand S. 1996. Cyclical behavior of thrust wedges: Insights from high basal friction sandbox experiments., 24(2): 135–138.

Gutscher M, Kukowski N, Malavieille J, Lallemand S. 1998. Material transfer in accretionary wedges from analysis of a systematic series of analog experiments., 20(4): 407–416.

Itasca Consulting Group. 2008. PFC2D (Particle Flow Code in 2 Dimensions) Online Manual (Version 4.0). Minnesota: Itasca Consulting Group Inc.

Jaeger J C, Cook N G, Zimmerman R. 2009. Fundamentals of Rock Mechanics. John Wiley & Sons.

Li C S, Yin H W, Jia D, Zhang J X, Wang W, Xu S J. 2018. Validation tests for discrete element codes using single- contact systems., 18(6), 6018011.

Li C S, Yin H W, Wu C, Zhang Y C, Zhang J X, Wu Z Y, Wang W, Jia D, Guan S W, Ren R. 2021. Calibration of the discrete element method and modelling of shortening experiments., 9, 636512.

Li C S, Zhang D, Du S S, Shi B. 2016. Computed tomography based numerical simulation for triaxial test of soil-rock mixture., 73: 179–188.

Liu Z, Koyi H A, Swantesson J O, Nilfouroushan F, Reshetyuk Y. 2013. Kinematics and 3-D internal deformation of granular slopes: Analogue models and natural landslides., 53: 27–42.

Malvern L E. 1969. Introduction to the Mechanics of a Continuous Medium. New Jersey, United States. Prentice- Hall, Inc.

Mathworks. 2015. MATLAB documentation. Natick, United States. The MathWorks Inc.

Means W D. 2012. Stress and strain: Basic concepts of continuum mechanics for geologists. Springer Science & Business Media.

Morgan J K. 2015. Effects of cohesion on the structural and mechanical evolution of fold and thrust belts and contractional wedges: Discrete element simulations.:, 120(5): 3870– 3896.

Oertel G. 1996. Stress and deformation: A handbook on tensors in geology. Oxford University Press.

O’Sullivan C. 2011. Particulate Discrete Element Modelling: A Geomechanics Perspective. Taylor & Francis.

Schreurs G, Buiter S J, Boutelier D, Corti G, Costa E, Cruden A R, Daniel J, Hoth S, Koyi H A, Kukowski N. 2006. Analogue benchmarks of shortening and extension experiments.,,, 253: 1–27.

Schreurs G, Buiter S J, Boutelier J, Burberry C, Callot J, Cavozzi C, Cerca M, Chen J, Cristallini E, Cruden A R. 2016. Benchmarking analogue models of brittle thrust wedges., 92: 116–139.

Sun C, Jia D, Yin H W, Chen Z X, Li Z G, Shen L, Wei D T, Li Y Q, Yan B, Wang M M, Fang S Z, Cui J. 2016. Sandbox modeling of evolving thrust wedges with different preexisting topographic relief: Implications for the Longmen Shan thrust belt, eastern Tibet.:, 121(6): 4591–4614.

Suppe J. 2007. Absolute fault and crustal strength from wedge tapers., 35(12): 1127–1130.

Wang X Y, Morgan J, Bangs N. 2021. Relationships among forearc structure, fault slip, and earthquake magnitude: Numerical simulations with applications to the Central Chilean Margin., e2021GL092521. doi:10.1002/essoar.10506057.2

Xu W Q, Yin H W, Jia D, Li C S, Wang W, Yang G X, He W H, Chen Z X, Ren R. 2021. Structural features and evolution of the northwestern Sichuan Basin: Insights from discrete numerical simulations., 9, 653395.

Yin H W, Zhang J, Meng L S, Liu Y P, Xu S. 2009. Discrete element modeling of the faulting in the sedimentary cover above an active salt diapir., 31(9): 989–995.

Quantitative Analysis and Simulation of Compressive Tectonics Based on Discrete Element Method

LI Changsheng1, 2, 3, 4, YIN Hongwei3*, XU Wenqiao3, WU Zhenyun1, 2, GUAN Shuwei4, JIA Dong3, REN Rong4

(1. State Key Laboratory of Nuclear Resources and Environment, Nanchang 330013, Jiangxi, China; 2. School of Earth Sciences, East China University of Technology, Nanchang 330013, Jiangxi, China; 3. School of Earth Sciences and Engineering, Nanjing University, Nanjing 210023, Jiangsu, China; 4. PetroChina Research Institute of Petroleum Exploration & Development, Beijing 100083, China)

With development of discrete element theory and computer technology, the discrete element method has been widely used in structural simulation at different scales. Compared with the traditional scaled analog (sandbox) experiments, the discrete element numerical simulation can precisely control the experimental boundary conditions and quantitatively analyze structural deformation process, which is helpful for deeply understanding the mechanical properties and deformation mechanism of stratum from the mesoscopic scale. In this paper, based on a typical discrete element simulation of compressive tectonics, we systematically expounded the tectonic deformation quantitative analysis method, highlighted the superiority of tectonic deformation in mesoscale, and obtained the following results: (1) In the typical numerical simulation of compressive tectonics, the forward spreading thrust imbricated faults are dominant, and the faults generate in sequence from the compression end to those far from the compression end. (2) Fractures are closely related to the formation of faults, the interval in which the fracture generation reaches its peak precedes the one in which the fault appears the most active. (3) In the early stage of fault formation, the faults have a small displacement. However, the number of new fractures is the largest in this stage. In the later stage of fault activity, the fault-slip extension amount increases linearly with the shrinkage of the model, and the number of new fractures is less. (4) Volumetric strain can represent fracture types (tension or compression), and deformation strain can distinguish forward and back thrust faults. (5) The average stress is positively correlated with topographic relief, and the maximum shear stress continues to accumulate under the new fault to be formed until the new fault is formed, then the shear stress begins to dissipate and spread forward, accumulating at the next new fault to be formed. The results of this study show that the discrete element method has great potential in the quantitative research of structural deformation, stress-strain, and fracture prediction.

discrete element method; structural deformation; compressive tectonics; stress and strain; quantitative analysis

P542; TE352

A

1001-1552(2022)04-0645-017

10.16539/j.ddgzyckx.2022.04.001

2020-12-06;

2021-06-16

國(guó)家自然科學(xué)基金項(xiàng)目(42102270、41972219、42172232、41927802、41572187、41602208)、東華理工大學(xué)博士啟動(dòng)基金項(xiàng)目(DHBK2019024、DHBK2019053)和中國(guó)石油天然氣股份有限公司項(xiàng)目(2018A-0101)聯(lián)合資助。

李長(zhǎng)圣(1989–), 男, 博士, 主要從事構(gòu)造及巖土體相關(guān)數(shù)值模擬的研究工作。E-mail: sheng0619@163.com

尹宏偉(1971–), 男, 教授, 主要從事構(gòu)造分析和構(gòu)造物理、數(shù)值模擬等領(lǐng)域的研究工作。E-mail: hwyin@nju.edu.cn