一种基于有限元的岩石圈长期变形数值计算方法

杨少华; 李忠海

doi:10.12090/j.issn.1006-6616.2018.24.06.079

一种基于有限元的岩石圈长期变形数值计算方法

doi: 10.12090/j.issn.1006-6616.2018.24.06.079

杨少华^1,,
李忠海^2, ,

1.
中国地质科学院地质研究所深地动力学重点实验室, 北京 100037
2.
中国科学院大学计算地球动力学重点实验室, 北京 100049

基金项目:

国家自然科学基金项目 41604080

国家自然科学基金项目 41622404

国家自然科学基金项目 41774108

国家自然科学基金项目 41590860

博士后面上基金项目 2016M601083

中石油前陆专项 2016B-0501

详细信息

作者简介:
杨少华(1987-), 男, 助理研究员, 主要从事计算地球动力学研究。E-mail:yangshaohua09@sina.com

通讯作者:
李忠海(1982-), 男, 教授, 主要从事计算地球动力学研究。E-mail:li.zhonghai@ucas.ac.cn

中图分类号: P313
计量
- 文章访问数: 329
- HTML全文浏览量: 142
- PDF下载量: 25
- 被引次数: 0
出版历程
- 收稿日期: 2018-08-13
- 修回日期: 2018-08-31
- 刊出日期: 2018-12-01

A NUMERICAL CALCULATION APPROACH BASED ON FEM FOR LONG-TERM DEFORMATION OF LITHOSPHERE

YANG Shaohua^1
,,
LI Zhonghai^{2
, ,}

1.
Institute of Geology, CAGS, Key Laboratory of Deep-Earth Dynamics, Beijing 100037, China
2.
Key Laboratory of Computational Geodynamics, University of Chinese Academy of Sciences, Beijng 100049, China

摘要

摘要: 有限单元法以其灵活性和精确性，成为固体地球科学中广为使用的数值方法。从短周期的地震活动到长周期的岩石圈变形、地幔对流，甚至行星演化，有限单元法几乎在固体地球科学的各个领域都占据着十分重要的位置。随着研究的深入，某些特定的地学问题给有限元计算带来挑战，尤其是岩石圈尺度大变形的数值计算，比如俯冲带的演化、剪切带中塑性流变导致的应力集中。基于显式有限元，尝试考虑粘弹塑性岩石圈大变形过程的数值计算。应用Marker-In-Cell（MIC）方法处理物质迁移。在描述基本原理和流程的基础上，对粘弹性变形、弹塑性变形、大变形过程及热传递过程等核心模块分别做了基准测试，而这四个模块是模拟岩石圈长期变形的关键。由测试结果和其他学者的（解析或数值）研究结果比对情况来看，受测试的核心模块基本达到了测试要求。可以预见，现有的基本算法可以满足研究岩石圈大变形的需要，进一步的具体研究工作将探讨这类问题。从科学问题层面讲，逐渐复杂的科学问题有利于数值模型的成熟。已达到基准测试的数值方法对下一步开展一些具体的地球动力学数值模拟研究有实际意义。
- 有限单元法 /
- 大变形 /
- 粘弹塑性材料 /
- 地球动力学
Abstract: Finite element method (FEM), which is very important in almost all fields of solid earth sciences, has been widely used in numerical experiments in solid earth sciences due to its flexibility and precision, ranging from short-term seismic activity to long-term lithospheric deformation, mantle convection and even planetary evolution. However, with the deepening of the research, some specific geological problems bring challenges to the finite element calculation, especially the numerical calculation of large-scale lithosphere deformation, such as the evolution of subduction zone and stress localization caused by plastic rheology in shear zone. Based on explicit FEM, an attempt is made to numerical calculation in the process of large deformation of lithosphere considering visco-elastic-plasticity in this study. Marker-In-Cell (MIC) approach was used in processing each material migration. On the basis of describing the basic principle and process, the core modules of visco-elastic deformation, elastic-plastic deformation, large deformation and heat transfer were tested. These four modules are the key to simulate the long-term deformation of the lithosphere. According to the comparison between the test results and the previous simulation results, the tested core modules basically meet the test requirements. It can be predicted that the existing basic algorithms can meet the needs of studying the large deformation of the lithosphere, and further specific research work will explore this kind of problem.
- finite element method /
- large deformation /
- visco-elastic-plastic materials /
- geodynamics

HTML全文

0. 引言

岩石断裂既是地壳构造变形和演化的重要产物，也是一种重要的构造作用与构造变形方式，它对油气运移和油气藏形成以及寻找金属矿床等很多方面都具有重要意义^[1]。

关于岩石断裂的机理和准则近百年来国内外众多专家和学者都做过很多探索。经典断裂力学理论中的一个根本假设是将岩石断裂轨迹视为直线型平面模型，然而现场实测和实验观测均表明，无论是在晶粒尺度还是在断层尺度上，岩石的断裂面都是非常不规则和粗糙的，难以用一个平直面来近似模拟^[2~3]，显示出了经典断裂力学理论的局限性。除此之外常用的断裂准则还包括：库仑断裂准则、水平直线型莫尔包络线理论、抛物线型莫尔包络线理论、格里菲斯断裂准则^[4]。以卡门实验为例，抛物线型莫尔包络线理论更符合实际情况^[5]。

关于地壳中岩石断裂的剖面形态，国内外学者都做过很多研究和探讨，总结出了常见的主要断裂样式：平直型断裂（层）、铲式断裂（层）和坡坪式断裂（层）^[6]。这些认识符合自然界实际情况，得到了普遍的认可，能定性地解释和预测一些现象，但对断裂产生机理的定量解释和断裂形态的定量表达却很少^[7~8]。本文将在抛物线型莫尔包络线理论基础上，定量地探讨上拱力构造背景下产生正断裂的形态、断裂极限应力及其与上拱力、地层物性之间的关系。这将对于精确预测解释断裂、认识盆地、认识上拱构造以及预测油气藏、分析油气运移等有重大意义。

1. 抛物线型莫尔包络线理论

当莫尔包络线以抛物线来近似表示（见图 1）时，其方程可写为^[5]：

图 1 抛物线型莫尔包络线^[5]

τ, σ—断裂面上的临界剪应力和临界正应力，MPa；σ₁, σ₃—发生断裂时最大和最小主应力，MPa；σ_m—平均应力（σ_m=(σ₁+σ₃)/2），MPa；σ_I—岩石在各向等值拉伸条件下的抗张断裂极限，MPa；

$\frac{\theta }{2}$

—岩石剪裂角，(°)

Figure 1. The Parabolic Mohr failure envelope

下载: 全尺寸图片幻灯片

${\tau ^2} = \frac{{{K^2}}}{{{\sigma _{\text{I}}}}}\left( {{\sigma _{\text{I}}} - \sigma } \right)$

(1)

式中：σ_I为岩石在各向等值拉伸条件下的抗张断裂极限，MPa；K为岩石粘结强度，MPa；τ、σ分别为断裂面上的临界剪应力和临界正应力，MPa。本文中规定，正值表示拉力，负值表示压力。

以主应力来表达，抛物线型莫尔包络线准则可写成^[5]：

${\left( {{\sigma _1} - {\sigma _3}} \right)^2} + 2\frac{{{K^2}}}{{{\sigma _{\text{I}}}}}\left( {{\sigma _1} + {\sigma _3}} \right) = {K^2}\left( {4 - \frac{{{K^2}}}{{{\sigma _{\text{I}}}}}} \right)$

(2)

式中：σ₁、σ₃分别为发生断裂时最大和最小主应力，MPa。

按照这一准则，剪裂角可表达成^[5]：

$\left. \begin{gathered} \tan \theta = 2\frac{{{\sigma _{\text{I}}}}}{K}{\left( {1 - \frac{{{K^2}}}{{2{\sigma _{\text{I}}}}} - \frac{{{\sigma _1} + {\sigma _3}}}{{2{\sigma _{\text{I}}}}}} \right)^{\frac{1}{2}}} \hfill \\ \theta = \arctan 2\frac{{{\sigma _{\text{I}}}}}{K}{\left( {1 - \frac{{{K^2}}}{{2{\sigma _{\text{I}}}}} - \frac{{{\sigma _1} + {\sigma _3}}}{{2{\sigma _{\text{I}}}}}} \right)^{\frac{1}{2}}} \hfill \\ \end{gathered} \right]$

(3)

随着σ₁与σ₃以及K和σ_I的变化，剪裂角可由0°变到45°。从（3）式可以看出，剪裂角不仅与物性参数（K，σ_I）有关，而且与各点极限应力状态有关。即使在同一种岩石中，由于各处极限应力状态不同，也可以形成不同的剪裂角^[5]。下面以此理论为基础讨论上拱力构造背景下正断裂剖面形态和极限应力状态。

2. 正断裂剖面形态和极限应力状态的理论分析

研究表明，地壳深部一般应力状态为：

${{\sigma _1} = {\sigma _2} = {\sigma _3} = - \rho gh}$

(4)

式中：ρ表示岩体的密度，kg/cm³；g为重力加速度，m/s²；h表示参考点的深度，m；σ₁、σ₂、σ₃分别表示最大主应力、中间主应力、最小主应力，MPa。

裂陷盆地中正断裂的形成常与盆地深部岩浆或地幔柱上升等有关^[9~12]。在上拱力构造背景下，岩石除受到一般应力的作用，还受到由上拱构造引起的附加水平拉力和上拱力共同作用。这种情况下，假设岩体受到的上拱力恒定，为σ_y2，附加水平拉力为σ_x2，埋深为h点处岩体竖直方向一般应力为σ_y1，水平方向一般应力为σ_x1，埋深为h点处岩体竖直方向总应力为σ_y，水平方向总应力为σ_x，则有：

${{\sigma _{y2}} = - \rho g{h_2}}$

(5)

${{\sigma _{y1}} = - \rho gh}$

(6)

${{\sigma _{x1}} = - \rho gh}$

(7)

${{\sigma _y} = {\sigma _{y1}} + {\sigma _{y2}}}$

(8)

${{\sigma _x} = {\sigma _{x1}} + {\sigma _{x2}}}$

(9)

式中：σ_y2为上拱力，MPa；h₂为上拱力以一般应力形式表示时对应的深度，m，由σ_y2导出。

当该点处岩体承受的竖直和水平方向总应力达到极限应力状态而发生正断裂时，根据抛物线型莫尔包络线理论，将（2）式变形整理可得：

${{\sigma _1} = {\sigma _3} - \frac{{{K^2}}}{{{\sigma _{\text{I}}}}} \pm 2K\sqrt {1-\frac{{{\sigma _3}}}{{{\sigma _{\text{I}}}}}} }$

(10)

在正断裂中，由于σ_x＞σ_y，所以有：

${{\sigma _1} = {\sigma _x},{\sigma _3} = {\sigma _y}}$

(11)

联立（5）—（11）式，整理可得深度为h点处岩体发生正断裂的极限应力状态：

${{\sigma _y} = - \rho g\left( {h + {h_2}} \right)}$

(12)

${{\sigma _x} = - \rho g\left( {h + {h_2}} \right) - \frac{{{K^2}}}{{{\sigma _{\text{I}}}}} + 2K\sqrt {1 - \frac{{ - \rho g\left( {h + {h_2}} \right)}}{{{\sigma _{\text{I}}}}}} }$

(13)

分析（13）式中水平方向极限应力并结合水平方向一般应力大小，σ_x由两部分组成，分别是地壳水平方向一般应力σ_x1和由上拱构造引起的附加水平拉力σ_x2，它们都随深度而变化。其中：

${{\sigma _{x2}} = - \rho g{h_2} - \frac{{{K^2}}}{{{\sigma _{\text{I}}}}} + 2K\sqrt {1 - \frac{{ - \rho g\left( {h + {h_2}} \right)}}{{{\sigma _{\text{I}}}}}} }$

(14)

当岩体所承受的总应力达到断裂极限应力状态发生断裂（见图 2）时：

图 2 岩体正断裂剖面数学分析

l—衡量地面水平长度的量，m；h—参考点深度，m；

$\frac{\theta }{2}$

—岩石剪裂角，(°)；σ_x, σ_y—深度为h点处岩体水平方向和竖直方向总应力，MPa

Figure 2. Mathematical analysis of geometric shape of normal fault of rock mass

下载: 全尺寸图片幻灯片

$\frac{{dl}}{{dh}} = \tan \frac{\theta }{2}$

(15)

式中l为衡量地面水平长度的量，m。

联立（3）、（11）、（12）、（13）、（15）式，可求得断裂剖面的数学表达式为：

$\frac{{dl}}{{dh}} = \frac{{\sqrt {1 + \frac{{{K^2}}}{{4\sigma _{\text{I}}^2}} + \frac{{\rho g\left( {h + {h_2}} \right)}}{{{\sigma _{\text{I}}}}} - \frac{K}{{{\sigma _{\text{I}}}}}\sqrt {1 + \frac{{\rho g\left( {h + {h_2}} \right)}}{{{\sigma _{\text{I}}}}}} } - \frac{K}{{2{\sigma _{\text{I}}}}}}}{{\sqrt {1 + \frac{{\rho g\left( {h + {h_2}} \right)}}{{{\sigma _{\text{I}}}}} - \frac{K}{{{\sigma _{\text{I}}}}}\sqrt {1+\frac{{\rho g\left( {h + {h_2}} \right)}}{{{\sigma _{\text{I}}}}}} } }}$

(16)

式（16）即为在上拱力背景下，岩体受到上拱力σ_y2和由上拱构造引起的附加水平拉力σ_x2而达到极限应力状态发生正断裂时，断裂剖面的数学表达式。

分析（16）式可知，断裂的剖面形态和极限应力状态只与地层物性参数K、σ_I、ρ有关；而断裂的分布深度则与上拱力大小相关。

3. 正断裂剖面形态及其分布深度和断裂极限应力状态模拟

3.1 正断裂剖面形态模拟

根据对式（16）的分析，下文将研究上拱力大小为0，即在水平拉张背景下，断裂形态和极限应力状态。

中国石油大学（北京）对兰聊断裂带附近砂岩进行测定及计算得到：K=50 MPa、σ_I=16 MPa、ρ=2.32 g/cm³。用此砂岩建立各向同性均质砂岩体模型，模拟水平拉张背景下正断裂剖面形态和极限应力状态。根据理论分析，并将具体数值代入（12）、（13）、（14）及（16）式，可得到在水平拉张背景下，发生断裂时断裂的极限应力状态为：

${{\sigma _y} = - 0.02273h}$

(17)

${{\sigma _x} = - 0.023h + \sqrt {14h + 10000} - 156.25}$

(18)

要达到断裂极限应力状态需要提供的附加水平拉力为：

${{\sigma _{x2}} = \sqrt {14h + 10000} - 156.25}$

(19)

断裂剖面函数表达式为：

$\frac{{dl}}{{dh}} = \frac{{\sqrt {0.0014h - 3.1 \times \sqrt {0.0014h + 1} + 3.44} - 1.56}}{{\sqrt {0.0014h - 3.1 \times \sqrt {0.0014h + 1} + 1} }}$

(20)

根据所求得的断裂剖面函数表达式（20）画出断裂剖面（见图 3），可以看出：断裂从h₀=6169 m深度点开始出现，随深度增加，断裂由陡变缓，其内摩擦角由90°减小，无限接近0°，即剪裂角由0°开始，无限接近45°方向，最终在地下不够提供断裂极限应力的某点消失，呈铲式，对应于抛物线型莫尔包络线中应力圆与抛物线相切于两点的情况；对于从地表到地下h₀=6169 m深度点之间，从函数角度来讲，该表达式无解，显示出了该表达式的缺陷。而从实际地质力学角度来讲，存在竖直形态的断裂，其剪裂角为0°，对应于抛物线型莫尔包络线中应力圆与抛物线相切于其顶点的情况；深度为h₀的点，则对应于应力圆与抛物线顶点相切的极限情况。在本文中，将深度为h₀这种应力圆刚好与抛物线顶点相切、断裂开始出现弯曲的深度命名为起弯深度，将这个深度点命名为起弯断点（见图 3）。

图 3 正断裂剖面及对应的抛物线型莫尔包络线模型

l—衡量地面水平长度的量，m；h₀—起弯深度，m

Figure 3. The geometric shape of normal fault and parabolic Mohr failure envelope models

下载: 全尺寸图片幻灯片

3.2 正断裂剖面形态和极限应力状态与地层物性关系模拟

根据对（16）式的分析，断裂的分布深度与上拱力大小有关，上拱力越大，断裂分布越浅。下文将在上拱力（σ_y2=-ρgh₂）为-140.26 MPa（此值可使下文均质砂岩中的断裂起弯深度为0，其计算方法后文将阐述）的前提下模拟不同物性地层产生的不同断裂形态和极限应力状态。

3.2.1 均质砂岩中正断裂的剖面形态和极限应力状态模拟

根据对兰聊断裂带附近砂岩的测定及计算结果（K=50 MPa、σ_I=16 MPa、ρ=2.32 g/cm³）建立各向同性均质砂岩体模型，模拟此情况下正断裂剖面形态和极限应力状态。

将具体数值代入公式（12）、（13）、（16），可得此情况下极限应力状态和断裂剖面形态表达式：

${{\sigma _y} = - 0.022736\left( {{h_2} + 6169} \right)}$

(21)

${{\sigma _x} = - 0.022736h + \sqrt {100 + 0.14 \times \left( {h + 6169} \right)} - 296.5}$

(22)

断裂剖面的函数表达式为：

$\frac{{dl}}{{dh}} = \frac{{\sqrt {0.0014h - \sqrt {0.14h + 94.1} + 12.08} - 1.56}}{{\sqrt {0.0014h - \sqrt {0.14h + 94.1} + 9.64} }}$

(23)

由所求得的断裂剖面函数表达式画出断裂剖面（见图 4a），可以看出：在此上拱力的背景下，断裂的起弯深度恰好为0，即断裂恰好从地表开始角度发生变化，随深度增加，断裂剪裂角由0°开始，无限接近45°方向，最终在地下不够提供其断裂极限应力的某点消失。这与实际情况中铲式断裂能够很好地吻合，能够很好定量解释其形成机制。

图 4 上拱力构造背景下岩石极限应力及正断裂剖面形态

Figure 4. Critical stress state of rock and geometric shape of fault under tectonic setting of upwelling force

下载: 全尺寸图片幻灯片

3.2.2 均质泥岩中正断裂的剖面形态和极限应力状态模拟

取兰聊断裂带附近泥岩，测定并计算其参数，建立均质泥岩体模型。应用同样的步骤可以计算出其极限应力状态和断裂剖面形态的函数表达式，并画出断裂剖面（见图 4b），其形态与上述均质砂岩中断裂类似，不同的是该剖面起弯深度为-152 m，表示断裂从地上152 m开始剪裂角发生变化，而实际情况是断裂剖面只为其地下部分。此外，泥岩断裂剖面比砂岩剖面剪裂角变化更快，显示为更平缓。

3.2.3 均质砂泥岩互层断裂剖面形态和极限应力状态模拟

在上述均质砂岩模型中插入一薄层均质泥岩，建立砂泥岩互层模型。此时，断裂的极限应力状态和断裂剖面为同等情况下均质砂岩和泥岩模型中应力状态和断裂的组合（见图 4c），可以看出：该砂泥岩互层模型中，砂岩段断裂面相对较陡，泥岩段断裂面较缓，由深至浅组合而形成典型的“陡—缓—陡”的坡坪式正断裂，能够很好地定量解释实际情况中的坡坪式断裂的形成机制和其形态。

3.2.4 实际地层断裂剖面形态和极限应力状态模拟

实际情况下，包括K、σ_I、ρ在内的各项地层参数都发生变化。资料显示，这3个参数都随地层深度增加而增大，并且增大的速度越来越慢。通过近似模拟可得出：在此上拱力背景下，实际地层的断裂剖面形态与均质岩体模型中的形态类似，呈铲式断裂。

在已知K、σ_I、ρ三个参数随深度变化的关系时，可定量得出实际地层断裂的极限应力状态和断裂剖面的函数表达式。

3.3 正断裂分布深度与上拱力关系的模拟与分析

从数学角度分析（16）式可知，断裂分布的深度或起弯深度与上拱力有一次线性关系。

由于起弯断点对应于应力圆恰好与抛物线型莫尔包络线相切于其顶点的情况（见图 3），且：

$\tan \theta > 0$

(24)

所以可联立（3）、（11）、（12）、（13）、（24）式，求解得出：

${h > \frac{1}{{\rho g}}\left( {{\sigma _{y2}} + \frac{{{K^2}}}{{{\sigma _{\text{I}}}}} - {\sigma _{\text{I}}}} \right)}$

(25)

进而可求得起弯深度h₀:

${{h_0} = \frac{1}{{\rho g}}\left( {{\sigma _{y2}} + \frac{{{K^2}}}{{{\sigma _{\text{I}}}}} - {\sigma _{\text{I}}}} \right)}$

(26)

式（26）即可表示起弯深度随上拱力的变化关系。

因此，对于在上拱力背景下产生的断裂，断裂形态只与地层物性有关，断裂分布的深度或起弯深度与上拱力呈一次线性关系，上拱力绝对值越大，断裂分布深度越浅，反之则越深。当上拱力为0时，即为水平拉张背景下产生的断裂（见图 3）；当上拱力达到σ_y2=σ_I-K²/σ_I时，断裂起弯深度为0，断裂刚好从地表开始发生弯曲（见图 4a、图 5）

图 5 起弯深度随上拱力变化关系

σ_y2—上拱力，MPa；h—深度，m

Figure 5. The relationship of fault's depth and upwelling force

下载: 全尺寸图片幻灯片

4. 结论

在上拱力背景下，假设上拱力恒定，均质地层受到一般应力、上拱力和由上拱构造引起的附加水平拉力，当上拱力和附加水平拉力达到一定值时，地层处于极限应力状态发生正断裂，断裂极限应力状态和剖面函数表达式可以定量得出。这将对于断裂解释预测、极限应力状态和上拱力以及附加水平拉力求取、油气藏预测及分析具有重大意义。

断裂由浅到深可以分为两个部分：从地表到起弯断点之间，断裂呈现竖直形态；起弯断点以下，随深度增加，断裂由陡变缓，剪裂角由0°无限趋近45°，一直延伸至不够提供其断裂所要求的极限应力的地下某点消失，呈现出铲式断裂的形态。这与实际情况中对铲式断裂的认识相一致，能够很好地定量解释其形成机制。

断裂的形态只与地层物性有关，泥岩中产生的断裂剪裂角比砂岩角度变化快；断裂分布的深度与上拱力呈一次线性关系，上拱力绝对值越大，断裂分布深度越浅。

当在较硬均质砂岩中夹入较软泥岩层时，可形成典型的“陡—缓—陡”的坡坪式正断裂，断裂在硬岩层中较陡，软岩层中较缓。断裂剖面形态的函数表达式和极限应力状态亦可定量表达出来。能够很好地解释和预测坡坪式断裂的形态及其形成机理。

图 1 直角坐标中的三角形单元(e)的面积坐标

Figure 1. Triangle' s area-coordinates in cartesian coordinates

下载: 全尺寸图片幻灯片

图 2 马克斯韦尔体应力积累过程的边界条件和数值解与解析解的对比

a—边界条件；b—数值解(线)与解析解(圈)的对比

Figure 2. Stress accumulation of a Maxwell body: boundary conditions and the comparison between numerical solutions (line) and analytical solutions (circle)

下载: 全尺寸图片幻灯片

图 3 弹塑性体应力积累过程的边界条件和数值解与解析解的对比

a—边界条件；b—数值解(线)与解析解(圈)的对比

Figure 3. Stress accumulation of elastic-plastic body: boundary conditions and the comparison between numerical solutions (line) and analytical solutions (circle)

下载: 全尺寸图片幻灯片

图 4 不同作者对瑞利—泰勒不稳定性的测试结果对比

a—文中结果；b—Thieulot结果^[27]；c—Van Keken结果^[26]

Figure 4. Rayleigh-Taylor instability verifications

下载: 全尺寸图片幻灯片

图 5 岩石圈一维热传导过程中温度剖面数值解(线)和解析解(十字)的对比

Figure 5. 1-D temperature profile comparisons between numerical solutions and analytical solutions of the lithosphere

下载: 全尺寸图片幻灯片

参考文献(30)

[1]	Hughes T J R. The finite element method:linear static and dynamic finite element analysis[M]. New York:Dover Publication, 2000.
[2]	Ismail-Zadeh A, Tackley P. Computational methods for geodynamics[M]. Cambridge:Cambridge University Press, 2010.
[3]	Tucker G E, Hancock G R. Modelling landscape evolution[J]. Earth Surface Processes and Landforms, 2010, 35(1):28~50. doi: 10.1002/esp.v35:1
[4]	Yang S H, Shi Y L. Three-dimensional numerical simulation of glacial trough forming process[J]. Science China Earth Sciences, 2015, 58(9):1656~1668. doi: 10.1007/s11430-015-5120-8
[5]	Beaumont C, Jamieson R A, Nguyen M H, et al. Himalayan tectonics explained by extrusion of a low-viscosity crustal channel coupled to focused surface denudation[J]. Nature, 2001, 414(6865):738~742. doi: 10.1038/414738a
[6]	Van Keken P E, Kiefer B, Peacock S M. High-resolution models of subduction zones:implications for mineral dehydration reactions and the transport of water into the deep mantle[J]. Geochemistry, Geophysics, Geosystems, 2002, 3(10):1~20. doi: 10.1029/2001GC000256/pdf
[7]	Leng W, Zhong S J. Viscous heating, adiabatic heating and energetic consistency in compressible mantle convection[J]. Geophysical Journal International, 2008, 173(2):693~702. doi: 10.1111/gji.2008.173.issue-2
[8]	杨少华, 许志琴, 李忠海, 等.欧罗巴星陨石坑对冰层厚度的制约[J].地球物理学报, 2017, 60(3):935~940. http://www.cnki.com.cn/Article/CJFDTotal-DQWX201703009.htm YANG Shaohua, XU Zhiqin, LI Zhonghai, et al. Constraint of impact craters on ice thickness on the Europa[J]. Chinese Journal of Geophysics, 2017, 60(3):935~940(in Chinese with English abstract) http://www.cnki.com.cn/Article/CJFDTotal-DQWX201703009.htm
[9]	Zhang S Q, O'Neill C. The early geodynamic evolution of Mars-type planets[J]. Icarus, 2016, 265:187~208. doi: 10.1016/j.icarus.2015.10.019
[10]	He Y, Puckett E G, Billen M I. A discontinuous Galerkin method with a bound preserving limiter for the advection of non-diffusive fields in solid Earth geodynamics[J]. Physics of the Earth and Planetary Interiors, 2017, 263:23~37. doi: 10.1016/j.pepi.2016.12.001
[11]	Enright D, Losasso F, Fedkiw R. A fast and accurate semi-Lagrangian particle level set method[J]. Computers & Structures, 2005, 83(6~7):479~490. http://www.ams.org/mathscinet-getitem?mr=2143508
[12]	Weinberg R F, Schmeling H. Polydiapirs:multiwavelength gravity structures[J]. Journal of Structural Geology, 1992, 14(4):425~436. doi: 10.1016/0191-8141(92)90103-4
[13]	Gerya T V. Introduction to numerical geodynamic modelling[M]. Cambridge, UK:Cambridge University Press, 2010.
[14]	Ranalli G. Rheology of the Earth[M]. UK:Chapman and Hall, 1995.
[15]	Kronbichler M, Heister T, Bangerth W. High accuracy mantle convection simulation through modern numerical methods[J]. Geophysical Journal International, 2012, 191(1):12~29. doi: 10.1111/gji.2012.191.issue-1
[16]	Kaus B J P. Factors that control the angle of shear bands in geodynamic numerical models of brittle deformation[J]. Tectonophysics, 2010, 484(1~4):36~47. doi: 10.1016/j.tecto.2009.08.042
[17]	Moresi L, Dufour F, Mühlhaus H B. Mantle convection modeling with viscoelastic/brittle lithosphere:numerical methodology and plate tectonic modeling[J]. Pure and Applied Geophysics, 2002, 159(10):2335~2356. doi: 10.1007/s00024-002-8738-3
[18]	Kaus B J P, Popov A A, Baumann T S, et al. Forward and inverse modelling of lithospheric deformation on geological timescales[A]. NIC Symposium 2016[C]. NIC Series, 2016, 298~307.
[19]	Gerya T V, Yuen D A. Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems[J]. Physics of the Earth and Planetary Interiors, 2007, 163(1~4):83~105. doi: 10.1016/j.pepi.2007.04.015
[20]	Jaeger J C, Cook N G W, Zimmerman R W. Fundamentals of rock mechanics[M]. Malden:Blackwell Publishingl, 2007.
[21]	王勖成.有限单元法[M].北京:清华大学出版社, 2003. WANG Xucheng. Finite element method[M]. Beijng:Tsinghua University Press, 2003. (in Chinese)
[22]	夏志皋.塑性力学[M].上海:同济大学出版社, 1991. XIA Zhigao, Plasticity mechnics[M]. Shanghai:Tongji University Press, 1991. (in Chinese)
[23]	Shewchuk J R. Delaunay refinement algorithms for triangular mesh generation[J]. Computational Geometry, 2002, 22(1~3):21~74. doi: 10.1016/S0925-7721(01)00047-5
[24]	Burden R L, Faires J D. Numerical analysis[M]. Boston:Brooks/Cole, 2011.
[25]	李庆扬, 王能超, 易大义.数值分析[M]. 5版.北京:清华大学出版社, 2008. LI Qingyang, WANG Nengchao, YI Dayi. Numerical analysis[M]. 5th ed. Beijing:Tsinghua University Press, 2008. (in Chinese)
[26]	Van Keken P E, King S D, Schmeling H, et al. A comparison of methods for the modeling of thermochemical convection[J]. Journal of Geophysical Research:Solid Earth, 1997, 102(B10):22477~22495. doi: 10.1029/97JB01353
[27]	Thieulot C. ELEFANT:a user-friendly multipurpose geodynamics code[J]. Solid Earth Discussion, 2014, 6(2):1949~2096. doi: 10.5194/sed-6-1949-2014
[28]	Turcotte D L, Schubert G. Geodynamics[M]. Cambridge, UK:Cambridge University Press, 2002.
[29]	le Pourhiet L, May D A, Huile L, et al. A genetic link between transform and hyper-extended margins[J]. Earth and Planetary Science Letters, 2017, 465:184~192. doi: 10.1016/j.epsl.2017.02.043
[30]	Ruh J B, Gerya T, Burg J P. 3D effects of strain vs. velocity weakening on deformation patterns in accretionary wedges[J]. Tectonophysics, 2014, 615~616:122~141. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=a43bd9c299bbc09c6d683626a1ea36c3

被引

资源附件(0)

访问统计

图(5)

计量

文章访问数: 329
HTML全文浏览量: 142
PDF下载量: 25
被引次数: 0

0. 引言
1. 抛物线型莫尔包络线理论
2. 正断裂剖面形态和极限应力状态的理论分析
3. 正断裂剖面形态及其分布深度和断裂极限应力状态模拟
3.1 正断裂剖面形态模拟
3.2 正断裂剖面形态和极限应力状态与地层物性关系模拟
3.3 正断裂分布深度与上拱力关系的模拟与分析
4. 结论

0. 引言
1. 抛物线型莫尔包络线理论
2. 正断裂剖面形态和极限应力状态的理论分析
3. 正断裂剖面形态及其分布深度和断裂极限应力状态模拟
3.1 正断裂剖面形态模拟
3.2 正断裂剖面形态和极限应力状态与地层物性关系模拟
3.3 正断裂分布深度与上拱力关系的模拟与分析
4. 结论

参考文献(30)

被引

资源附件(0)

访问统计

姓名
邮箱
手机号码
标题
留言内容
验证码

留言板

一种基于有限元的岩石圈长期变形数值计算方法

doi: 10.12090/j.issn.1006-6616.2018.24.06.079

作者简介:
杨少华(1987-), 男, 助理研究员, 主要从事计算地球动力学研究。E-mail:yangshaohua09@sina.com

通讯作者:
李忠海(1982-), 男, 教授, 主要从事计算地球动力学研究。E-mail:li.zhonghai@ucas.ac.cn

计量