基于離散元的擠壓構造定量分析與模擬

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

基于離散元的擠壓構造定量分析與模擬

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

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

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

離散元; 構造變形; 擠壓構造; 應力應變; 定量分析

0 引 言

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

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

離散元數值模擬方法在構造形態、斷裂預測及應力應變的定量分析方面具有明顯的優勢。眾多學者已經將離散元引入構造變形的模擬中, 但分析的重點多集中于構造幾何形態的定性解析。本文在闡述離散元原理、模擬過程的基礎上, 結合一個典型的擠壓構造模擬實例, 詳細剖析了基于離散元的擠壓構造實驗過程及定量分析方法。

1 基本原理

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

相互粘結的兩個顆粒, 其接觸力的計算可以分為兩種狀態, 即壓縮與拉伸。當兩個顆粒處于壓縮狀態下, 接觸力計算采用Hertz-Mindlin接觸力學模型; 當兩個顆粒處于拉伸狀態下, 接觸力計算采用粘結模型。

一般地, 規定兩個顆粒間壓力為正, 拉力為負。

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

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

(1) 壓縮狀態

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

其中,

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

當顆粒處于壓縮狀態下時,K由下式計算,

(2) 拉伸狀態

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

圖1 接觸面定義

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

當顆粒處于粘結狀態時, 粘結的抗拉力為,

粘結的抗剪力為,

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

當兩個顆粒處于非粘結狀態(斷裂)時,Fmax= 0.0, 顆粒間不產生拉力;Fmax=μF, 切向力需滿足|F|≤Fmax。

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

2 應力應變表征

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

2.1 應 變

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

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

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

體積應變I計算方法如下:

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

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

2.2 應 力

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

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

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

3 實驗過程

3.1 參數標定

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

3.2 擠壓實驗

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

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

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

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

4 結果分析

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

圖2 雙軸實驗結果

表1 顆粒的細觀參數

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

表2 顆粒間的粘結參數

注：壓縮速度取2.0 m/s, 實驗過程是準靜態過程, 詳細調試過程及結果見李長圣(2019)。

圖3 初始模型

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

4.1 楔體構造定量分析

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

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

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

4.2 斷裂解譯

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

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

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

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

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

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

4.3 應力和應變

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

圖7 每擠壓1 km過程中, 新斷開的粘結力鏈數

收縮量不同時, 模型內的累積變形應變見圖11。變形應變強的區域與體應變較高的區域相互重疊(蔡申陽等, 2016)。正向逆沖(紅色)往往伴隨有反向逆沖(藍色)斷層, 尤其當底部滑脫層強度較弱時, 這種共軛的“X”狀發育的斷層更為明顯(Morgan, 2015;李長圣和尹宏偉, 2017; 李長圣, 2019)。斷層F3(圖11c)和斷層F4(圖11d)的淺部變形應變較大, 也表明斷層是從淺部開始發育, 并逐步向深部發展。

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

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

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

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

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

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

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

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

5 結 論

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

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

(2) 裂縫與斷層形成有密切關系, 局部區域內聚集的大量裂縫是產生斷層的誘因;

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

圖14 收縮量為20 km時,水平方向平均應力分布

黑色線條為應力等值線, 黑色點為|變形應變|>4的點, 表征了斷層的位置。圖(e)中, 綠色層為標志層。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

吳珍云, 尹宏偉, 李長圣, 孫業君, 李麗梅, 杜航, 何奕成, 劉博雅. 2019. 斷陷盆地正反轉構造實驗模擬及對茅山斷裂帶的啟示. 南京大學學報(自然科學版), 55(5): 869–878.

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

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

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

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

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

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

李長圣(1989–), 男, 博士, 主要從事構造及巖土體相關數值模擬的研究工作。E-mail: sheng0619@163.com

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