李长圣, 尹宏伟, 徐雯峤, 吴珍云, 管树巍, 贾 东, 任 荣

基于离散元的挤压构造定量分析与模拟

李长圣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, 接着楔体继续增厚, 坡角显著增大, 进入下一个演化循环。新断层开始活动, 楔体高度持续增大, 楔体宽度相对稳定, 坡角减小。应力累计、断裂持续增加, 形成断层, 使得坡角发生波动, 这种波动值大小可能与地层强度有关。总体上, 坡角趋于一个稳定值, 符合临界楔理论(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