构造应力场模拟——有限元理论、方法和研究进展①

张胜利

（1.五邑大学信息学院，广东江门 529020；2.中科院广州地化所，广东广州 510640）

构造应力场模拟
——有限元理论、方法和研究进展①

张胜利1，2

（1.五邑大学信息学院，广东江门 529020；2.中科院广州地化所，广东广州 510640）

采用有限元数值模拟方法对构造地质问题进行描述和定量化求解是当前地质学领域的研究的一个热点，在近10年以来取得了重要进展，形成了比较完整的理论和技术体系，并在一些典型的地质构造带获得了重要的研究成果。本文以有限元数值模拟方法理论作为出发点，总结分析了国内外有限元数值模拟方法在构造应力场领域的研究进展情况和技术方法，并讨论了其目前存在的问题和未来发展方向。

构造应力场；数值模拟；有限单元法

Abstract：The Finite Element Method（FEM）has been used in the study of tectonic stress field for a long time，and the essence of numerical modeling has been adopted to the well－established numerical methods of multidisciplinary acknowledge including mathematics，physics and mechanics for studing characters of geological tectonics.In the last decade，great advances have been made on the numerical simulation method，not only an integrated theory has been built up，but also some significant results have been born from several typical tectonic belts.So the FEM becomes one of the most important numerical methods in the study of tectonic stress field.In this paper，taking theory of FEM as a springboard，the new progress and methods in this field at home and abroad is summarized and analyzed.Some problems and prospect of the researching on the field is also given.

Key words：Tectonic stress field；Numerical model；Finite element method

0 引言

地壳中的各种地质构造都是岩石受力发生变形的产物，它们的产生和发展必然也受力学规律的支配。因此研究构造应力场特征是阐明和解释各种地质构造的产生、分布和演化规律的重要途径。但由于地质构造演化过程的特殊性，采用传统的地质－地球物理－地球化学的方法研究应力场特征显得异常困难。而物理模拟虽然可以帮助人们理解构造变形过程，但是存在着严重的时空尺度的局限性，并且一些因素在实验室条件下尚不能获得，尤其是在伴随着高温高压的地球深部构造带。

有限单元法又称有限元法，是上世纪50年代中期发展起来的一种数值计算方法，最初主要用于结构和应力分析，60年代后期被引入到地质构造分析中，由于其巨大优越性而被迅速推广，已成功的解决了很多地学的难题。我国是从70年代中期开始将有限元数值模拟技术应用在地学领域的。随着有限元数值模拟技术的出现，在综合利用其它方法（如地质、地球物理、地球化学）研究成果的基础之上，采用这类数值模拟技术，可以建立不受时空限制的各种地质模型，帮助研究人员比较直观地认识和理解构造应力场的特征和演化历史，并可以对其演化趋势做出预测。现在该技术除了在基础地质方面得到广泛应用外，在天文、环境、矿产勘探、地震预测等方面也取得了众多的研究成果［1－3］。

随着计算机计算速度的提高，构造应力场的数值模拟方法也获得了很大的进展。国内外的研究结果表明，目前已经可以在构造应力场的模拟过程，在不同的本构关系下，将流体、温度、剥蚀、断层等因素考虑进去，进行二维和三维的数值模拟［4－5］。但是也存在很多的问题，例如模拟的主观性（边界条件、模型等）、模拟结果的多解性和可靠性等。

1 基本理论

有限元法的基本原理是根据三大守恒定律（质量、动量、能量）建立三个偏微分方程：平衡方程、本构方程和几何方程。然后将拟分析的连续体假想地分为有限个单元，离散后的单元与单元之间利用单元节点相互连接起来。单元内部点的待求量可由单元节点量通过选定的插值函数求得。由于单元形状简单，易于用平衡关系或能量关系建立节点量之间的方程式。然后将各个单元方程“组装”在一起形成总体代数方程组，计入边界条件后即可对方程组求解。单元划分越细，计算结果就越精确。

采用有限元法模拟构造应力场一般包括一下几个步骤：

（1）在详细分析研究目标地质资料的基础之上建立计算模型，这也是数值模拟试验成败的关键。

（2）单元剖分和选择插值函数。根据模型的几何特性、载荷情况及所要求的变形点，将模型离散化，单元大小应当可以给出有用的结果，但又要尽可能的节省计算量。

（3）确定应变位移和应力应变关系。应力与应变之间关系由本构方程确定，不同的本构关系对应不同的物理介质。目前构造应力场模拟中较常用的本构关系有4种：线弹性、弹塑性、粘弹性和弹粘塑性本构关系［6］。

其中弹塑性本构一般用于模拟地壳浅部的脆性层，粘弹性和弹粘塑性本构多用于模拟地壳深部。也可以在同一个模型中采用多种本构关系。

（1）线弹性本构关系：

其中σij是应力张量；eij是应变张量；λ和μ是拉梅系数；θ是体积模量。线弹性本构一般用于模拟小范围的地壳浅部沉积层的应力场。

（2）弹塑性本构关系：

其中Dijkl为弹性系数；εkl为应变张量；dλ是由一致性条件确定的尺度因子；F是用应变空间表示的屈服函数。弹塑性本构一般用于模拟地壳浅部脆性层。

（3）粘弹性本构关系（以Kelvin体为例）：其中E为弹性模量；εij为应变速率张量；η是粘性系数。式（3）也常被称为线性流变方程。

（4）弹粘塑性本构关系：

其中γ是控制塑性流动速率的流度参数，是时间和粘塑性应变张量的不变量的函数；f是屈服函数；f0是一个使表达式无量纲化的、正的参考值。公式（2）、（3）也可以视作本式的特殊形式。

2 研究历史及现状

构造应力场模拟和板块碰撞、大陆变形数值模拟紧密相关。Bott［7］对板缘力驱动下板内应力应变的演化模拟可以说是这方面的开山之作。1976年Tapponnier［8］等人用理想刚塑性体所产生的变形来模拟亚洲大陆自新生代以来所经历的大尺度碰撞变形和断层的产生，并为板块内部强烈变形提供了有力证据，从而使得应力场的数值模拟方法引起了人们的广泛关注。England［9］等在1982年首次将粘性薄层流变模型运用到以印度板块和欧亚板块的碰撞为代表的大陆变形带，并考虑了地壳增厚造成的应力变化，对青藏高原的挤压隆升演化进行了数值模拟；England和Houseman等［10－11］也成功模拟了印度板块和欧亚板块碰撞后板块边缘变形与地幔对流之问的耦合关系，并指出地壳厚度变化是造成应力场变化的主要因素，走滑断裂只是部分地缓解了这种应力场的变化。

但是上述研究中也存在明显的不足之处，如考虑的只是小变形，忽略了温度、重力、沉积、剥蚀等因素的影响等。但实际的地质过程中，一方面地质构造过程属于大变形（几何非线性），并不满足上述有限元法变形（小变形）的前提条件，此外所有的地质因素是相互影响耦合在一起的，单独考虑某一个和几个因素必然会影响模拟的结果，甚至出现大的偏差。因此Houseman和England［10－11］对大变形，大旋转的速度场和应力场等进行了较全面的分析，并提出了新的有限元模拟方法。Zhang等用有限元的模拟手段成功地解释了澳大利亚东部的岩石圈构造和异常应力场的成因［12］，并用变形一流体一热力学一化学反应全耦合的方法实现了在构造控矿理论研究和矿产勘探中的应用。

实际上以加拿大人Beaumont为核心的研究小组已经建立了一个完整的数值模拟体系，包括理论方法，计算技术以及各类计算模［13－15］等。他们采用不同的边界条件进行计算，对斜向汇聚碰撞和巨型造山带的应力场进行了模拟分析，模拟介质包括弹塑性，刚塑性，粘弹性和粘弹塑性等，同时还考虑了温度场、应力场和流体场的耦合问题，并较好的解释了新英格兰阿巴拉契亚、喜马拉雅、阿尔卑斯等著名造山带的应力异常问题［5，13－16］。此外，Beaumont C、Batt G E等人在研究挤压构造带时，还全面的分析了山体抬升引起的附加重力、剥蚀过程对质量的重新分配，以及由此造成的应力应变的变化。并给出了构造变形与剥蚀过程的耦合数值求解方法，其中计算剥蚀量的地表剥蚀模型包括了长距离的河流搬运作用和短距离的岩崩、滑坡、冲刷等扩散作用的网格地形模型［17－18］。

国内一些地质学者在有限元数值模拟方面也做了很多重要工作。其中青藏高原由于其在全球构造中独特性，以及有限元数值模拟方法在研究构造应力场演化机制方面的优越性，使得青藏高原一直是国内专家学者的研究热点和重点，同时也是世界上数值模拟方法应用最多的一个区域之一。

石耀林［19］等用有限元方法研究分析青藏高原热演化中长期变形粘性流动，是国际上最早对陆陆碰撞非弹性力学分析的论文之一。熊熊、傅容珊等［20－21］假设地幔是粘滞度为常数的二维牛顿流体，模拟了青藏高原增厚大陆岩石层热边界被对流地幔剥离并为软流圈物质替代的动力学过程，并推断高原地壳和岩石层在岩石层增厚和剥离以前就很热，指出挤压隆升受多种因素的影响，而且高原隆升是不均匀演化的过程。郑勇等［22］采用有限元方法模拟了近20万年来青藏高原岩石圈形变演化过程，探讨了印度一欧亚大陆的碰撞对中国大陆岩石层形变和应力场的影响以及它们与强地震活动性的关系。段岩等［23］人结合深部地质与地球物理资料，对喜马拉雅构造带中的札达盆地断裂的构造应力场进行了模拟，计算结果表明札达盆地的演化明显受盆地两侧边界断裂的控制，札达盆地是在整体南北向挤压应力的作用下不同块体差异隆升作用的结果。曹建玲等［24］根据GPS测量结果，利用三维球坐标系下的Maxwell体有限元模型模拟了青藏高原在印度板块推挤下的变形。他们认为柔软下地壳的存在使整个青藏高原在印度板块的推挤下表现为整体抬升，高原南缘和喜马拉雅东构造结抬升迅速，而周边地块下地壳相对较硬而封闭，仅高原东南部存在高温柔软的通道，使得高原整体隆起达到一定高度后，下地壳和软流层的物质向东、东南流动，并拖曳上地壳作类似运动，形成绕喜马拉雅东构造结的顺时针转动。

此外，有限元数值模拟方法在地震成因以及孕震过程等方面也得到了非常广泛的应用，这方面的研究已有许多报道。早在1972年，Lysmer［25］就将有限元法应用于地震数值模拟之中，Mat suura等［26］用一个简化了的岩石圈－软流圈耦合模型数值模拟了横穿板块边界的构造载荷情况，研究结果表明地震断层的应力累积一部分来源于其基底载荷（基底岩石圈的粘性阻力），一部分来源于其边缘载荷（断层水平边缘上的位错累积）；而Hashimoto等［27－28］通过综合考虑粘弹性滑动响应函数、随时间与滑动变化的断层本构关系，以及给定一个板块运动驱动力建立了仿真3D模型，计算模拟了板块边界地震发生的整个过程。在国内，王仁等人［29］是较早开展这方面的研究人员。朱守彪［30］以2004年发生过特大地震的苏门答腊俯冲带为例，模拟了俯冲带上俯冲板片与上伏板块之间的闭锁、解锁、滑动到再闭锁这一准周期性过程，成功的模拟了地震的孕育、发生过程。实际上目前有限元数值模拟方法孕震过程与前兆机理研究中的应用已经取得了很大进展：地质模型从二维到三维；本构关系从单一的弹（塑）性到粘－弹性、粘－塑性，从线性到非线性；模型中地壳（岩石圈）的分层从简单的上、下两层到细划分多层［1－2，31－34］；计算方法从原来的单机串行计算到目前的CPU矩阵并行计算甚至是云计算等［35］（更多的文献可在有限元方法在地球科学中应用的评述性文章中查到［36］）。刘启元、马宗晋［36］等人曾经指出，地震预报研究应充分吸收和借鉴气象预报的成功经验，以GPS为代表的空间对地观测技术，岩石圈巨型高分辨率宽频带流动地震台阵观测技术以及数值模拟技术已经为地震数值预报研究提供了前所未有的技术基础。当前地震数值预报突破的目标无疑应首先主要集中在地震的中期或中短期预报。以地震数值预报为目标的GPS阵列地壳形变连续观测、高分辨率地壳上地幔结构探测、地壳动力学、地震孕育和破裂过程的理论、模拟试验和实际观测、数据同化和计算软件的开发应成为今后研究的重点。地震数值预报研究必将极大地促进中国地震科学基础研究和地震预报的进一步发展。

尽管国内的研究人员已将构造应力场的有限元数值模拟技术广泛地应用在了地震、地质、油田开发等诸多行业，并且在计算方法、模型处理做了很多有益的尝试，但是在模型建立、技术处理等方面与国际上还存在着一定的差距，并且目前还没有开发出能够获得同行普遍认同的地学专业有限元软件。很多研究者还是采用的线性模型，用一些通用的有限元软件（ANSYS，ALGOR等）来处理复杂地质问题，存在较大的误差，致使一些模拟结果难以令人信服，此外也较少考虑时间以及地表过程（沉积、剥蚀、滑坡等因素）的影响等。

3 相关的技术方法

由于数值模拟结果存在多解性，因而建立一个合适的计算模型是模拟构造应力场中十分重要的一个环节，从而减少多解性。它包括建立地质模型和数值（物理）模型两个阶段：在建立数值模型阶段，需要定义模型的几何形状，选择表达模型物理行为特征的数学描述方法，并确定适合该地质模型有限单元类型；在建立地质模型阶段，需要根据观测或推测的边界条件、实验参数，来建立实验研究的初始条件和边界条件，并在模拟过程中不断地修正、更改实验参数、初始条件和边界条件进行模拟实验，逐步使实验结果趋近于观察到或科学推测的地质构造模型，从数学理论和方法上来证实或证伪已有的地质模型，从而建立起自己的符合区域地质实况的地质模式。因此，综合分析已有资料，构建合理的地质模型，选择不同的数学模型进行实验研究，并从不同方面选择实验研究的初始条件、边界条件及实验参数，从地质、地球物理、地球化学等方面提供约束条件，限制实验结果的多解性，并使实验结果趋近于地质现状，是数值模拟的主要研究内容与方法。

此外，在构造应力场的模拟过程中，如果考虑地壳深部温度场的影响，还需要引入热动力学方程，将热动力学方程和三大方程（平衡方程、本构方程和几何方程）进行联合求解，即温度场和应力场的耦合处理。Fullsack［13］、Beaumont［15］、Batt［18］等对热一力学耦合问题给出了相应的结果。

采用有限元数值方法模拟构造应力场的过程中，断层的处理是其中至关重要的一环。因为有限元法较适合分析连续介质的力学问题，对于地质构造中涉及到的节理、断层等不连续现象需要特殊的处理方法。一般断层问题可以采用以下几种办法解决：一种是goodman单元法［37］，这种单元比较简单，但是不适合算大的位移（变形）；一种是滑移节点法，这种方法常用在模拟地震的瞬间演化；另外是拉格朗日乘子法，这种方法可以直接求出断层面上的接触力，因为里面涉及一个满矩阵的求逆，所以计算规模比较小，但它的最致命的弱点是不容易收敛。对于地质上大的构造变形问题（如造山运动），如果本构关系选择密律流体本构，断层的处理会相对容易些，只需要把断层单元划分薄一点，然后给断层单元分配一个比较低的粘度系数就可以了。并且选择密律流体本构也优于粘弹性本构，因为大变形都涉及长时间的应力松弛现象。这也是现在比较通行一种处理方法［17］，熊熊、曹建玲等［20－24］均采用了该方法处理构造应力场中的断层。

当要考虑地表过程造成的影响时，例如由于山体抬升引起的附加重力、剥蚀沉积过程对质量的重新分配，以及由此造成的应变应力的变化等，需要采用特殊的耦合处理方法，有的还需要在剖分网格的时候进行局部加密。Beaumont［15－16］等给出了构造变形与剥蚀过程的耦合数值求解方法。但是该耦合过程不包含基本变量的耦合求解，只是在每一个时间步的构造抬升计算完成之后，根据剥蚀模型计算剥蚀量，再对地形进行修正。其中计算地壳变形抬升的模型为受挤压作用的刚塑性或粘性平面应变模型；计算剥蚀量的地表剥蚀模型为包括了长距离的河流搬运作用和短距离的岩崩、滑坡、冲刷等扩散作用的网格地形模型。

4 问题与讨论

用有限元法模拟构造应力场不仅是当前定量研究地质力学的热点和焦点，也是难点之一。虽然我国在这方面的研究工作早已经开始，但是力量薄弱，只能说还处于探索阶段，与国外同行还有不小的差距，尤其是多场耦合的三维非线性问题的模拟。此外构造应力场的模拟涉及的学科领域较多，需要考虑固体力学、流体力学和热力学等各方面的信息，对于沉积、剥蚀等地表因素的处理也都需要用不同的方程表示，并需要将上述各方面的控制方程耦合在一起，这样就极大的增加了求解难度。要很好的解决这些问题需要大的学科交叉和科学家的合作，尤其需要数学家、物理学家的参与，并尽早的开发出适用于地学的有限元专业软件。

［1］ 章纯.中国东部地区地震活动与构造应力场关系的有限元数值模拟［J］.西北地震学报，2007，29（3）：230－234.

［2］ 王爱国，袁道阳，梁明剑.兰州盆地最大潜在地震变形数值模拟［J］.西北地震学报，2008，30（3）：232－238.

［3］ 王建，叶正仁.地幔对流对全球岩石圈应力产生与分布的作用［J］.地球物理学报，2005，48（3）：584－590.

［4］ 陆远忠，叶金铎，蒋淳，等.中国强震前兆地震活动图像机理的三维数值模拟研究［J］.地球物理学报，2007，50（2）：499－ 508.

［5］ Beaumont C，Jamieson R A.Himalayan tectonics explained by extrusion of a low－viscosity crustal channel coupled to focused surface denudation［J］.Nature，2001，414：738－742.

［6］ 殷有泉.计算固体力学方法［M］.北京.科学出版社，2003：121－148.

［7］ Bott M H，Dean D S.Stress systems at young continental margins［J］.Nature，1972，235：23－25.

［8］ Tappennier P，Molnar P.Slip line field theory and large－scale continental tectonics［J］.Nature，1976，264：319－324.

［9］ England P C，Mackenzie D A.Thin viscous sheet model for continental deformation［J］.Journal of Geophysics Research，1982，70：295－321.

［10］ England P，Houseman G.Finite strain calculation of continental deformation，2.Comparison with the India－Asia collision zone［J］.Journal of Geophysical Research，1986，91（B3）：3664－3676.

［11］ England P，Houseman G.Role of Lithospheric strength hetemgeneities in the tectonics of Tibet and neighboring regions［J］.Nature，1985，315：297－301.

［12］ Zhang Y，Scheibner E，OrdA.Numerical modeling of crustal stresses in the eastern Australian passive margin［J］.Austra1ian Journal of Earth Sciences，1996，43：161－175.

［13］ Fullsack P.An arbitrary Lagrangian－Eulerian formulation for creeping flows and its application in tectonic models［J］.Geophysical Journal International，1995，120：1－23.

［14］ Ellis S，Beaumont C.Models of convergent boundary tectonics：implications for the interpretation of Lithoprobe data［J］.Earth Sciences，1999，36：1711－1741.

［15］ Beaumont C，Ellis S，Hamilton J.Mechanical model for subduction－collision tectonics of Alpine－type compressional orogens［J］.Geology，1996，24：675－678.

［16］ Beaumont C，Muoz J A，Hamilton，et al..Factors controlling the Alpine evolution of the central Pyrenees inferred from a comparison of observations and geodynamical models［J］.Journal of Geophysical Research，2000，105：8121－8145.

［17］ Willett S，Beaumont C，Fullsack P.Mechanical model for the tectonics of doubly vergent compressional orogens［J］.Geology，1993，21：371－374.

［18］ Batt G E，Braun J.On the thermo－mechanical evolution of compressional orogens［J］.Geophysical Journal International，1997，128：364－382.

［19］ Chi－yue Wang，Yao－lin Shi.Dynamic uplift of the Himalaya［J］.Nature，1982，298：553－556.

［20］ 熊熊，傅容珊，许厚泽，等.增厚大陆岩石圈热边界层对流剥离的数值模拟［J］.地球物理学报，1998，41（增刊）：33－40.

［21］ 傅容珊，徐耀民，黄建华，等.青藏高原挤压隆升过程的数值模拟［J］.地球物理学报，2000，43（3）：346－355.

［22］ 郑勇，傅容珊，熊熊.中国大陆及周边地区现代岩石圈演化动力学模拟［J］.地球物理学报，2006，49（2）：4153－427.

［23］ 段岩，孟宪刚，邵兆刚，等.西藏札达盆地控盆断裂有限元数值模拟［J］.地质通报，2008，27（10）：12－16.

［24］ 曹建玲，石耀霖，张怀，等.青藏高原GPS位移绕喜马拉雅东构造结顺时针旋转成因的数值模拟［J］.科学通报，，2009，54（2）：224－234.

［25］ Lysmer J.A Finite element method for seismology［J］.Journal of Computational Physics，1972，（11）：181－216.

［26］ Mat suura M，Sato T.Loading mechanism and scaling relation of large interplate earthquakes［J］.Tectonophysics，1997，277：189－198.

［27］ Hashimoto C，Mat suura M.3Dsimulation of earthquake generation cycles and evolution of fault constitutive properties［J］.Pure Appl Geophys，2002，159：2175－2199.

［28］ Aochi H，Mat suura M.Slip and time2dependent fault constitutive law and it s significance in earthquake generation cycles［J］.Pure.Appl.Geophys.，2002，159：2 02922 046.

［29］ 王仁，孙荀英，蔡永恩.华北地区近年地震序列的数学模拟［J］.中国科学（B辑），1982，18（8）：715－753.

［30］ 朱守彪，邢会林，谢富仁，等.地震发生过程的有限单元法模拟—以苏门答腊俯冲带上的大地震为例［J］.地球物理学报，2008，51（2）：460－468.

［31］ 王爱国，石玉成，柳煜.黄河黑山峡大柳树坝址区最大潜在地震变形及地震应力模拟预测［J］.西北地震学报，2007，29（4）：314－318.

［32］ 文彦君.延怀盆地设定地震浅析［J］.西北地震学报，2008，30（2）：159－162

［33］ 刘海明，陶夏新，孙晓丹，等.马衔山北缘断裂西段6.5级地震对兰州市及周边地区的影响［J］.西北地震学报，2008，30（3）：227－231.

［34］ 何丽君，石玉成，杨惠林，等.地震动作用下黄土边坡稳定性分析［J］.西北地震学报，2009，31（2）：142－147.

［35］ 张怀，吴忠良，张东宁，等.虚拟川滇—基于千万网格并行有限元计算的区域强震演化过程数值模型设计和构建［J］.中国科学（D辑），2009，39（3）：260－270.

［36］ 刘启元，吴建春.论地震数值预报—关于我国地震预报研究发展战略的思考［J］.地学前缘，2003，10（8）：217－224.

［37］ Goodman R E.Methods of geological engineering in discontinuous rocks［M］.Paul：West Publishing Co.，1976.

Modeling of Tectonic Stress Field——the Theroy，Method and Related Research Progress of the Finite Element Method

ZHANG Sheng－li1，2
（1.School of Information，Wuyi University，Guangdong Jiangmen 529020，China；2.Guangzhou Institute of Geochemistry，Chinese Academy of Sciences，Guangzhou 510640，China）

P315.12

A

1000－0844（2010）04－0405－06

2009－07－20

中国科学院知识创新工程项目（KZCXZ－SW－117）；广东高校优秀青年创新人才培养计划项目（LYM10128）

张胜利（1979－），男（汉族），东单县人，博士，主要从事数值模拟及数据挖掘方面研究.