李長圣, 尹宏偉, 徐雯嶠, 吳珍云, 管樹巍, 賈 東, 任 榮

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

李長圣1, 2, 3, 4, 尹宏偉3*, 徐雯嶠3, 吳珍云1, 2, 管樹巍4, 賈 東3, 任 榮4

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

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

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

0 引 言

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

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

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

1 基本原理

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

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

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

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

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

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

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

其中,

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

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

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

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

圖1 接觸面定義

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

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

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

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

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

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

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

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

2.1 應(yīng) 變

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

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

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

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

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

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

2.2 應(yīng) 力

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

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

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

3 實驗過程

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

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

3.2 擠壓實驗

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

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

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

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

4 結(jié)果分析

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

圖2 雙軸實驗結(jié)果

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

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

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

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

圖3 初始模型

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

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

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

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

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

4.2 斷裂解譯

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

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

黑色線條為應(yīng)力等值線, 黑色點為|變形應(yīng)變|>4的點, 表征了斷層的位置。圖(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)力開始消散, 繼續(xù)往前傳播, 在下一個將要形成的新斷層處累積。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

鄒瑋, 余一欣, 劉金水, 蔣一鳴, 唐賢君, 陳石, 余浪. 2021. 東海盆地西湖凹陷中央反轉(zhuǎn)構(gòu)造帶發(fā)育主控因素及寧波背斜形成過程. 石油學(xué)報, 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

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

李長圣(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