断层系统摩擦动力学行为的有限元模拟分析

邢会林，郭志伟，王建超，张熔鑫，刘骏标，姚琪

1 深海圈层与地球系统前沿科学中心,海底科学与探测技术教育部重点实验室, 中国海洋大学海洋高等研究院/海洋地球科学学院,青岛 266100 2 青岛海洋科学与技术国家实验室 海洋矿产资源评价与探测技术功能实验室,青岛 266100 3 中国海洋大学海底科学与工程计算国际中心,青岛 266100 4 中国地震局地震预测研究所，北京 100036

0 引言

地震预报是地球科学界的世界性难题，通过基于地球物理场观测、岩石物理实验和定量化的物理规律的数值模拟，研究地震孕育、破裂的物理机理，指导地震危险性的有效评估，最终实现基于物理的地震数值预报已成为人们的共识(石耀霖等,2013,2018;Barbot et al.,2012;Segall,2012;刘启元和吴建春，2003).王仁等(1980，1982)曾尝试利用有限元方法模拟分析华北地区地震序列及迁移规律，随后朱岳清等(1988)和蔡永恩等(1999)曾分别就唐山地震孕震和震源动力过程数值分析进行了初步的有益尝试.石耀霖等(2018)在总结统一的美国加州地震破裂预测系统基础上，提出了我国地震数值预报路线图的设想，吹响我国地震数值预报的起床号.

世界上80%以上的大地震发生在板块边界及板块内部断层上，被认为是断层或板块边界上的一种摩擦失稳行为(Brace and Byerlee,1966).中国大陆具有以板内块体为单元的构造变形特征，且强震多发生在活动块体的边界/断层上.地震的孕育和发生过程也就是断层闭锁、解锁和滑动过程，其中断层面上的接触摩擦失稳行为起着决定性的作用，为此国内外学者对此进行了大量的研究.一般而言，含断层的有限元分析主要有两种不同的处理方式：一种把断层当作某种特殊材料(如弱材料)的连续介质(如王仁等,1982;陈化然等,2004;王凯英和马瑾,2004;Zhang et al.,2009;Luo and Liu,2010;He et al.,2011;Liu et al.,2016;Li et al.,2017;Sun and Liu,2018)，另一种把断层作为非连续的摩擦接触界面(如Goodman,1968;Melosh and Williams,1989;Huang et al.,1997;Xing and Makinouchi,2002a,b;朱守彪和张培震,2009;Xing and Xu,2011;姚琪等,2012a,b,2018a;吴萍萍,2014;范桃园等,2014).两者都可计算相关的应力及变形行为，但后者在理论上或可更加真实地模拟断层的动力学行为.基于摩擦接触的断层有限元分析，摩擦面法向上多采用罚函数、拉格朗日乘子法或其扩展方法处理.‘Slippery node’方法(Melosh and Williams,1989)和双节点节理模型(Goodman,1968)一样，实质上采用了类似拉格朗日乘子法的处理都引入了新的未知变量，且其特殊处理也仅适用于二维模型计算(Huang et al.,1997).目前含断层的有限元分析多集中在长时间变形分析或库仑应力场计算，断层的摩擦系数多取为常数，摩擦失稳细节基本不考虑，且多为二维单一断层模型(Huang et al.,1997;朱守彪和张培震,2009;范桃园等,2014;吴萍萍等,2014).然而在现实中，大量的地震观测和研究表明地震孕育、破裂并非单一断层的孤立行为，而是以系统的形式存在.一方面，地震可能是由若干断层级联式破裂的，例如2016年新西兰地震破裂了10条左右断层.另一方面，同震破裂、震后应力调整使得周围断层面应力增加或减弱，激发大量余震或延迟触发大地震.断层系统的有限元数值模拟是了解它们的动力学和复杂的系统行为的有力工具，但是以往的研究由于观测资料少、断层系统精细网格难以生成等原因导致所用的模型过于简单粗糙，加之断层系统数值模拟技术不够稳定和成熟等客观条件的限制，对地震过程的数值模拟研究还十分有限.

虽然也有一些内部和商业软件用于模拟一些与地震有关的摩擦现象，但含断层相互作用的地壳动力学系统仍然是一个非常复杂的问题，因为沿断层的黏着-滑移摩擦失稳具有很强的非线性.目前，绝大多数基于摩擦接触的断层/地震行为的研究多采用机械工程行业的商业化软件，如ANSYS(如吴萍萍等,2014)和ABAQUS(如范桃园等,2014)来模拟断层行为，且摩擦系数多简单地取为常数.但是，机械工程行业常采用润滑剂降低界面摩擦系数，避免加工面发生划痕损伤同时降低能耗，摩擦系数一般最大为0.1左右.而地震断层上的静态摩擦系数一般在0.6以上，而且摩擦系数在摩擦失稳/地震发生时会急速下降，比如到0.2，而非一个常数.由此导致的高度非线性问题使得目前的商业化软件(如ANSYS、ABAQUS、MSC-MARC和ADINA等)计算稳定性遇到了巨大的挑战.此外，有的研究者尝试用动态显式法来模拟地震的破裂过程，但是计算时间步长太小，对于地震的孕育过程及整个地震周期的模拟则需要太长的计算时间，应力计算精度也多有问题.表1参照较为成熟的典型变形过程有限元计算分析，比较了不同时间积分方案的异同，其规律对地震过程的模拟也基本适用.

表1 静态隐式、静态显式和动态显式的比较Table 1 Comparison of static implicit,static explicit and dynamic explicit formulations

在目前的有限元地震孕育过程模拟分析中，静态隐式算法曾被尝试应用于此类非线性接触问题.表2对此进行了部分总结.有限元软件静态隐式算法通常采用Newton、Newton-Raphson或者修正的Newton-Raphson迭代法，其将控制方程作为一个非线性方程(组)，计算相应的位移增量来求解，迭代到平衡为止.然而Newton-Raphson等迭代算法的一个主要缺点是迭代可能会不收敛，因为对于非线性问题每步开始前难以给出足够精确的初始解.更何况，为了模拟地震行为，还要处理断层上的状态变化(如接触与脱离、黏着与滑移)以及滑动速率、滑移及接触压力等相关的摩擦系数变化等新的突变性非线性问题，而传统商业化软件面临这些问题往往收敛性极差.虽然采用迭代平滑技术来减少剧烈变化，但收敛仍然无法得到保证，而且这可能会导致黏着-滑移失稳发生时刻(即地震发生时刻)不能被准确地捕捉到.

表2 准静态摩擦接触问题有限元软件Table 2 Survey of some FEM software for quasi-static frictional contact problems

综上所述，虽然通过实验室试验和现场观测地震断层摩擦本构关系的研究取得了一系列成果，提出了诸如滑动位移、速度-状态相关等不同的非线性摩擦本构模型(Dieterich，1978，1979；Kato and Tullis,2001)，但是目前鲜有商业有限元程序可用于研究与地震孕育成核和破裂相关的非线性摩擦失稳、滑移动力学行为.现行代码/软件大多在接触界面上使用了预设的摩擦系数定值，而这些摩擦系数值与可变条件无关，因此无法研究变形岩石对摩擦界面行为的影响，如黏着-滑移摩擦失稳和沿接触界面的非线性滑动摩擦行为.虽然或许可以通过用户子程序来实现地震过程模拟所需的部分功能，但是由于如上所述非线性问题的收敛性难题可能导致不收敛而无解；如果采用迭代平滑技术来减少剧烈变化从而满足收敛条件，会导致黏着-滑移失稳(即地震发生的时刻)等关键突变时刻漏失而不能被准确捕捉到.

因此，Xing等(1998)、Xing和Makinouchi(2002a，2002c)提出了一种基于R-minmium自适应控制时间步长的静态显式时间积分方案，兼顾了静态隐式与动态显式的部分优点(见表1和表2中PANDAS)，可以精确捕捉黏着-滑移失稳等关键时刻；利用节点-点接触单元，将其应用于可变形体之间的三维任意形状非线性摩擦接触问题的有限元模拟，进一步利用速率和状态相关的摩擦定律分析断层的非线性摩擦失稳行为，取得了较好的效果，具体总结如下.

1 摩擦接触有限元算法

1.1 有限元公式

为了更好地描述地震过程尤其是局部大滑移大变形等问题，我们采用更新的拉格朗日公式来描述有限变形问题.结合物体间摩擦接触条件，其平衡方程可以采用虚拟速度原理等效表示为(Xing et al.,1998，Xing and Makinouchi,2002a,b,c)：

(1)

(2)

(3)

当材料表现为率相关，也就是黏弹性或黏弹塑性时，材料的本构方程可以表述为

(4)

1.2 摩擦接触的统一化本构方程

1.2.1 法向接触力

采用罚函数方法处理法向接触不可穿入的约束关系.当两个点在接触面的法线方向上的相对位移小于零时，即gn<0时发生接触，其法向接触力为：

fn=f·n=Engn(≠0仅限于gn<0),

(5)

式中，En为法向上的罚因子，gn为法向穿透量，f为接触力，n为接触体表面的外法线方向.

1.2.2 摩擦力

断层面上的摩擦接触关系是断层行为变化的关键因素，断层面的形状、运动状态、摩擦系数及其不均匀分布和演化等条件都可以通过断层的摩擦接触行为来表现，这为使用非线性摩擦有限元方法来进行区域地震危险性分析提供了可行性.

大多数构造地震都是由沿先存断层或板块界面的突然滑动引起的.大地震被认为是沿着断层的黏着-滑移摩擦失稳的结果(Brace and Byerlee，1966).此后，基于实验结果提出了几种摩擦本构方程，目前非线性摩擦问题可归纳为以下两种模式(Dieterich,1978，1979；Kato and Tullis,2001;Xing et al.,2006)，一种是滑动位移依赖型，一种是速度-状态依赖型.这两种模式都可以用来描述断层摩擦过程中的某种现象，但是关注点不一样.滑动依赖模式着重于破裂过程的描述，譬如从黏着状态转换到滑动状态的动态过程，但忽略了速度效应；而速度-状态依赖模式描述了滑动过程中状态以及滑动速度的影响，但难以描述从黏着/闭锁状态到产生初始滑动/破裂的行为.其实从实验结果也可以看出，两种模式描述了断层摩擦行为的不同阶段.在模拟断层行为和地震过程中，人们一直在试图将两种模式结合起来使用(Xing et al.,2006,2007).Xing等(Xing and Makinouchi,2002a,2002c,2003;Xing et al.,2006,2007)导出了统一化的数学公式表述速度相关的摩擦接触中黏着和滑移不同的运动状态，并导入有限元分析软件PANDAS，计算在内外载荷作用下接触面上摩擦接触状态(黏着和滑移两种状态)及变化过程，即断层上闭锁和解锁滑移两种不同运动状态的变化特征，也就是地震及其周期的孕育、发生的过程(参见表2中PANDAS).

(6)

(黏着状态,l=1,2)

(7a)

(7b)

1.2.3 速度-状态摩擦本构关系

Dieterich(1978，1979)和Ruina(1983)提出速率和状态相关摩擦关系式来描述断层的摩擦行为，已得到广泛应用；据此，与滑动速率相关的瞬态摩擦行为可表示为：

μ=μ0+φ+aln(V/Vref),

dφ/dt=-[(V/L)(φ+bln(V/Vref))].

(8)

把式(8)代入式(7)，也可以得到摩擦的类似材料本构式(4)的统一化表达式.此外，摩擦力是一个矢量，它也会随着接触面法线方向的变化而变化(Xing and Makinouchi,2002a,c;Xing et al.,2006).

1.3 时间积分方法

如前所述，为了应对摩擦接触状态、摩擦非线性以及材料非线性所引起的收敛性问题，我们采用了显式时间积分方案.假设在足够小的时间步长内，只要状态不发生剧烈变化(积分点发生弹塑性变化、摩擦界面节点发生接触脱离或者黏着-滑移状态变化)，公式(1)中的所有变化率在从t到t+Δt的时间增量内都可以被认为是恒定的.我们采用R-Minimum方法来限制步长，以避免在一个增量步长时间内状态发生剧烈变化，从而精确捕捉到状态变化的时刻(Xing and Makinouchi,2002a,2002c).

因此，式(1)中的所有变化率皆可以用增量来代替:

(9)

最后，结合上述公式(5)(7)和(8)，公式(1)有限元离散后可写为:

(K+Kf)Δu=ΔF+ΔFf,

(10)

式中，K和ΔF分别是整个物体的标准刚度矩阵和外载荷边界增量；当采用率相关材料模型时，它们也包括了黏性项的贡献(如，Xing and Zhang,2009;Xing and Wang,1997)；Kf和ΔFf是所有接触单元的接触刚度矩阵和接触力增量，Δu是节点位移增量(Xing and Makinouchi,2002a,c).

2 断层系统有限元计算网格生成

网格生成是进行有限元计算分析之前的关键一步，它是将研究对象离散为计算用网格(单元)，以便施加边界条件进行有限元模拟计算.目前主要使用商业化图形网格软件(如HyperMesh和Patran)以及ABAQUS、ANSYS等有限元分析软件中的前处理等功能模块进行网格生成.通过这些软件，可以将二维几何模型离散为三角形或四边形网格，或者将三维几何模型离散为四面体或六面体网格.由于这些商业化软件都是为机械和土木工程行业设计分析而开发的，一般通过确定的线、面和体等几何信息生成网格.这样的方法对于处理简单形状的地质模型尚可，但对于含复杂断层系统的地质体，一般在地质勘探图像中提取的相关断层数据并用一系列点阵来表示，而非确定的线、面和体等传统网格生成软件所需要的几何描述，因此含复杂断层系统地质体的计算用有限元网格模型生成还面临着较大的挑战.

目前常用的计算网格生成算法包括以下四种：Delaunay算法(Baker,1989;Borouchaki and Lo,1995;Borouchaki et al.,2000)，八叉树(四叉树)法(Frey et al.,1994;Schneiders,1996)，推进波前法(AFT)(Löhner,1996a,b;Lee and Hobbs,1999)以及扫描法或映射网格法(表3).为了在三维空间中生成包含由大量点数据集定义的断层网格，可以使用Delaunay算法或推进波前法自动生成四面体网格.但传统的Delaunay算法在应对复杂断层系统的有限元网格生成时并不理想，因为该方法并不能保证网格节点精确地插入边界或断层面处.也就是说，在生成的地质体网格中，断层面的形状和位置会发生变化，不能保证其原有的形状及特征.一般说来，虽然四边形和六面体网格比三角形和四面体网格具有更好的品质(Schneiders and Bünten,1995)，更适合于非线性有限元分析，但是对于三维地质体模型，特别是对于包含复杂断层的大规模地质体模型，目前仍难以自动生成高质量的非连续六面体网格(Lee and Lee,2002;Xing et al.,2009).

表3 表面域网格划分的主要算法及介绍Table 3 Main algorithms and introduction of surface mesh generation

2.1 基于断层面的非连续四面体网格生成

当断层系统的断层面非常复杂而曲折时，用六面体网格来对地质体模型进行网格划分不仅费时费力，而且往往只能将曲折界面简化为一系列具有简单几何形状的平面.为此，采用基于断层面约束的四面体网格自动生成方法，可以保证其中断层的几何形状、位置和整体网格的形状和质量，并用于有限元计算.具体的操作过程如下：首先对断裂带数据进行边界提取与约束，通过Delaunay与AFT算法相结合的方法建立空间三角形来描述断层空间曲面；其次对面网格进行优化处理，尤其是对断层曲面的交叉部分进行特殊识别及优化；接着为研究区生成一个由三角形网格组成的封闭外壳以完全包裹研究区断层面；最后以生成的内部断层及外壳的三角形面网格为约束条件，生成四面体网格(Liu and Xing,2016).这种方法的优点是可快速地自动生成含复杂断层系统的四面体有限元网格模型.然而，全自动化网格生成带来的缺点是生成的四面体网格单元和节点数量很多，网格质量尤其是含空间曲面断层时较难控制；同时对标准的线性四面体单元而言有限元应力计算精度较低，还在进一步优化研究之中.

2.2 非连续六面体网格生成

扫描法是非连续六面体网格生成较为常用的一种方法，通过直接或者非直接方法生成全四边形表面网格，然后以此按照预定的方向推进生成六面体形状的网格.同时将断层线推出的内表面拆分成两个单元间的分界面，形成沿断层的不连续摩擦接触界面.这种方法可以简单、快速地生成六面体网格，但生成的断层表面倾角多为90°，比较适用于走滑断层.这种方法的难点在于作为基础面的全四边形表面网格的生成上(Xing et al.,2009).Liu和Xing(2013)将三角形合成和波前推进法相结合，自动生成具有复杂线约束的高质量全四边形网格.下面，以含复杂裂隙结构材料的数字图像(图1a)为例，来说明利用该方法实现断层系统非连续六面体网格的生成.具体步骤如下：(1)断层结构提取与修复：利用PANDAS软件将图1a中的断层特征提取出来，并采用蚁群算法、DBSCAN聚类算法和Canny算子的混合图像特征修复算法(张愉玲等，2021)，恢复图像中缺失的断层特征(图1b)；(2)边界平滑：由于数字图像的性质，提取的边界几何形状通常是锯齿状的，并影响生成的四边形网格的质量，采用基于相位的边界平滑方法来平滑锯齿状边界，并对其进行离散，考虑到三维有限元模型进行应力加载和边界条件设置的边界效应，在断层区的四周加了一个较大的矩形边界，用于约束网格的生成和防止计算时边界出现突变(图1c).为了使断层周围的网格更加规则和精细，可以围绕断层生成一条或者几条尽量平行于断层的线作为新的约束线(图1d)；(3)初始三角形网格生成：利用改进的二维Delaunay三角剖分法生成具有断层线段约束的三角形网格(图1e)；(4)四边形网格生成：使用Catmull-Clark细分在每个三角形的质心处插入一个节点，将每个三角形分成三个四边形(将三角形网格转换为四边形网格).然后从线约束开始，基于生成的全四边形网格，通过推进波前法生成新的四边形单元；(5)四边形网格优化：通过优化四边形网格的拓扑结构以减少不规则节点的数量，并平滑生成的网格，交替执行全四边形网格的生成，直到生成高质量的网格(图1f)；(6)三维有限元网格生成：以四边形表面网格为基础种子面，沿着用户约定的方向推进扫描生成六面体形状的网格，然后采用潜在接触对方案对摩擦接触面进行拆解分离和标识(图1g).

图1基于图像生成的三维有限元模型(d,e,f中红色线段表示断层)(a)原始图像；(b)提取和修复后的断层；(c)平滑处理；(d)沿橙色断层生成蓝色等值线；(e)初始三角形网格生成；(f)四边形网格生成；(g)六面体网格生成(橙色四边形单元表示断层面).Fig.1 3D finite element model generated based on rock images (red line segments in d,e and f represent faults)(a)The original image;(b)Extracted and repaired faults;(c)Smoothing;(d)Blue contour lines generated along the fault in orange;(e)Initial triangular mesh generation;(f)Quadrilateral mesh generation;(g)Hexahedral mesh generation (Orange quadrilateral grids indicate fault planes).

对于断层在垂直方向上并不一定是竖直等特殊平面的情况，难以采用上述扫描法生成六面体网格，一般可采用映射分块法划分断层系统模型.本文以图2a所示的澳大利亚南部断层系统三维几何模型为例进行说明，其中部分断层面倾角非90°且并不是直线/平面.为了将该三维断层系统生成六面体网格并设定有限元模拟计算所必需的条件(即边界条件和关于断层的信息)，首先，在断层曲线/面的约束下，将整个断层系统几何模型分解成多个简单的可以网格化的六面体几何块体，其中断层面倾角为90°的可以用上述扫描法进行网格划分，其余的利用映射分块法进行网格划分(如图2b)；这些分组块各自生成六面体有限元网格后，在相邻组块体之间除断层接触面以外的公共边或曲面处共享相同的网格单元；最后，用“焊接”的方式(即删除重复点)或网格划分后用黏着接触算法将上述六面体网格组装在一起(如图3)，同时保持跨断层的网格不连续，在进一步的有限元分析中将其作为摩擦接触界面处理(Xing and Mora,2006;Xing et al.,2009).这种方法适用于可以将断层系统用映射分块法划分的模型，但对于复杂断层系统网格生成费时费力，且容易出错.一旦能够生成合理的网格，即可进行有限元计算分析.

图2 南澳断层系统(a)三维断层系统几何模型；(b)—(d)断层系统的不同几何组件/块体及其六面体网格;(e)组装在一起后生成的总体网格(修改自Xing and Mora,2006).Fig.2 The Southern Australian fault system(a)Geometry model of 3D fault system;(b)—(d)The hexahedron mesh generated from the different components/blocks of the fault system model;(e)The total mesh generated after assembling all the components together (Modified from Xing and Mora,2006).

3 摩擦接触有限元在断层系统分析中的应用

在上述的基础上，Xing等提出并采用R-Minimum自适应地控制不同非线性情况下的时间步长方案以保证计算精度和效率，将不稳定的隐式迭代计算转变为显式计算，可以自适应的控制时间步长以稳定高效地模拟地震孕育(时间可以为几十年甚至更长时间)和动态破裂(时间可能小到以分秒计)等地震过程及地震周期，采用基于R-Minimum的静态-显式方案开发了工业强度的摩擦接触有限元软件PANDAS(参见表2)，彻底解决了非线性断层系统摩擦接触的收敛性难题 (Xing and Makinouchi 2002a,b,c,2003;Xing et al.,2004,2006,2007).

目前，该软件已经成功模拟重现了经典的岩石摩擦实验，如三明治断层模型(Xing and Makinouchi,2002b)、单曲(Xing and Makinouchi,2003;Xing et al.,2004)和多曲的非平直三维断层模型(Xing et al.,2006).计算所得的接触面上的摩擦力、应变和应变率的演化与实验室发现的滑动弱化摩擦现象取得了很好的匹配.

针对现实中的多曲单断层摩擦行为模拟，该软件已经得到了很多的应用.Xing和Makinouchi(2002b),Xing等(2008)将这个方法应用到太平洋板块在日本海深部向下俯冲的模拟中，实现了连续曲度变化的三维断层面摩擦行为的收敛计算.Xing等(ACES会议报告，2006 5thACES,2008 6thACES,see http:∥www.aces.org.au)和朱守彪等(2008)利用该软件计算了多曲度苏门答腊俯冲带的摩擦行为，模拟了2004年苏门答腊地震的非线性孕震过程和地震破裂现象.姚琪等(2012a,b)利用该软件模拟计算了汶川地震起始破裂处的断层倾角和物性参数对强震的控制作用，证实了高倾角和上下盘巨大的岩性差异都能造成逆冲断层超长的孕震时间，并且龙日坝断裂分担了巴颜喀拉块体内部大部分的走滑分量，使得汶川地震起始破裂处以逆冲活动为主.姚琪等(2018b)利用该软件模拟了具有多曲度的喜马拉雅主逆冲带的断层破裂行为，用横向分段来解释了2015年尼泊尔地震复杂的地震破裂传播特征.郭婷婷等(2015)采用PANDAS对单断层与交叉断层两种断层模型分别进行数值模拟计算对比，结合中国大陆双震或震群型地震孕育发生的构造条件对共轭断层系统模型的孕震与发震机理进行了讨论与分析.

基于上述对非连续网格划分的方法和有限元基本理论的研究，本文以图1a中的复杂结构断裂系统为例进行研究.采用2.2节生成的非连续六面体网格模型作为本次模拟计算的复杂断层系统模型，三维有限元模型共有69595个节点，49356个八节点六面体(见图1f).这里使用速率相关摩擦定律(公式(8))来描述沿断层的摩擦行为，式中μ0=0.6，a=0.05，b=0.025;本文选取Vref作为整个过程的参考速度,同时假设模型中的材料具有相同的性质：密度ρ=2.60 g·cm-3，杨氏模量E=44.8 GPa，泊松比γ=0.22.模拟分析采用以下边界条件：在模型上面，沿Y轴的负方向以恒定的速度Vref/10施压，但Z方向位移固定；在底部，沿X、Y和Z三个方向固定；其他部分自由无边界条件约束.整个模型使用六面体单元进行三维计算，初始应力设为0.图3展示了不同阶段的X方向速度场(a—d)和总速度场(e—h)在XY平面分布变化情况.在加载开始阶段，数值模拟结果显示速度场及其X方向速度分量图表现为连续分布，断层两侧的速度分布显示还未发生明显的滑动(图3a和3e)；随着加载的进行，图3b和3f中与主加载方向(Y轴方向)夹角较小的断层，由于其断层面上切向摩擦力较大而法向正应力分量较小从而较易滑动，因此在较短的时间内由逐步闭锁状态进入局部滑动状态,并开始表现出局部断层区速度变化(图4a).随着持续加载，开始出现较大的断层破裂活动区，此时破裂区仍主要分布在低角度的断层部分(图3c和3g，图4b)；与主加载方向(Y轴方向)夹角较大的断层，在整体变形及周边断层滑动的影响下，也会逐渐进入滑动甚至脱离状态，导致断层两侧发生较大的速度突变(见图3d和3h,图4c)，断层系统的破裂范围发生进一步的扩大，整个计算区域整体明显向X轴的正负方向挤出.

图3 X方向速度场(a—d)和总速度场(e—f)在不同加载阶段的分布(单位：m·s-1).加载时间为100R/Vref其中(a)(e)R=0.001，(b)(f)R=0.003，(c)(g)R=0.014，(d)(h)R=0.06；深蓝曲线代表断层，(d)中黑色箭头表示块体挤出方向.Fig.3 The velocity component in X direction (a—d)and the total velocity magnitude (e—f)at different loading stages (unit:m·s-1).Loading time is 100R/Vref for (a)(e)R=0.001,(b)(f)R=0.003,(c)(g)R=0.014,(d)(h)R=0.06Dark blue curves represent faults,and black arrows in (d)represent the movement direction of blocks.

图4 (a)—(c)分别是图3(b)—(d)对应时刻间X方向速度分量的变化率深蓝色曲线代表断层.Fig.4 (a)—(c)are variation of velocity component in the X direction respectively between different moments in Fig.3 (b)—(d)Dark blue curves represent faults.

此外，该软件已经在地震相关的复杂断裂系统摩擦行为模拟中取得了多方应用，譬如Xing和Mora(2006)模拟了澳洲南部的复杂断层系统中(图2)断层的相互作用，证明了该软件能够用于不规整、多曲面、多断层、多接触的三维非线性摩擦模拟，实现了计算的收敛，获取应力场、应变场和断层节点摩擦力的时空演化特征，从而以此推测复杂断层系统中，多断层段逐次破裂的时空进程，图5显示了其中等效应力变化的一个特写镜头(Xing and Mora,2006).类似的，Xing等(2007)将美国南加州地区800 km×700 km×10 km范围内的9条断层进行了简化、建模和非线性摩擦有限元计算，模拟获取的等效应力场演化显示了不同规模断层的破裂、愈合、触发、延迟的复杂过程.针对国内川滇地区龙门山断裂、鲜水河—安宁河—小江断裂系、大凉山断裂等具有强震危险性的断裂，Xing和Xu(2011)、姚琪等(2018a)模拟了这些具有重要控制作用的边界断裂的破裂演化特征，模拟结果显示了在现今的应力场作用下，鲜水河—安宁河—小江断裂系和大凉山断裂的反复破裂-愈合过程，以及龙门山断裂带在长时间的应力加载下发生大规模破裂的过程.模拟所得的川滇地区主要断裂的破裂演化特征与该地区强震的时空分布基本一致，表明该软件所计算模拟的断层摩擦破裂行为能够代表强震的孕震-发震-震后愈合过程.

图5 变形某时刻等效应力变化率的分布快照(Xing and Mora,2006)Fig.5 A snapshot of the equivalent stress rate distribution at a certain deformation stage (Xing and Mora,2006)

4 结论与展望

地震预报是世界性难题，其重要价值不言而喻.要取得重大突破，基于物理的地震数值预报不可或缺.大地震被认为是断层上的一种摩擦破裂失稳行为，地震孕育和发生过程是断层系统相互作用动力学行为的一种表现，因此断层破裂摩擦动力学行为的研究对地震数值预报具有重要意义.本文围绕断层破裂摩擦动力学行为的有限元模拟，(1)概述了断层摩擦动力学有限元模拟相关的理论、软件及其优缺点，总结并比较了其中不同的时间积分方案和摩擦接触的计算方法，指出传统方案所面临的摩擦非线性的收敛性难题，即使采用迭代平滑技术来减少剧烈变化满足了收敛条件，也会导致黏着-滑移失稳(即地震发生的时刻)等关键剧变时刻漏失而不能被捕捉到；(2)提出并采用R-Minimum自适应地控制时间步长方案，将不稳定的隐式迭代计算转变为显式计算以保证各种非线性情况下的计算稳定性、精度和效率；同时通过节点-点接触单元，结合库仑摩擦准则利用速率和状态相关的摩擦系数分析断层非线性摩擦失稳行为；通过对相关应用分析总结，所提方案可自适应地控制时间步长以稳定高效地模拟地震过程，取得了较好的效果；(3)总结了基于断层系统有限元网格生成的不同方案，讨论了主流方法在复杂断层系统四边形及六面体网格生成中的弊端；通过复杂图像实例，介绍了所研发的基于图像的断层系统四边形和六面体网格自动生成步骤以及生成的网格模型在断层系统非线性有限元模拟计算中的应用.

为了提高断层系统有限元模拟分析平台的有效性和实用性，还需在以下几个方面进一步研究完善：

(1)多场耦合摩擦模型及其有限元实现：断层摩擦非线性行为十分复杂，除了受到接触面岩石物性、运动状态、滑动速度和位移、接触力的影响，还受到温度、流体孔隙压力和物理化学反应等各种耦合因素的影响.结合实验室试验和现场观测，发展断层摩擦多场耦合非线性本构方程及其稳定高效的有限元计算动力学模型与典型实验验证对于进一步探明多场耦合条件下摩擦失稳机理十分必要.

(2)多场耦合地震动力学综合研究：断裂带尤其是俯冲带是地球上构造活动最复杂、最强烈的区域.大洋岩石圈通过在汇聚板块边界的俯冲将大量水带入到地幔中，并对俯冲带地震的发生起到了重要的控制作用.随着俯冲板片深度的增加，在一定的温压条件下，水化地幔(即蛇纹岩)发生脱水相变，引发俯冲带中源地震.脱出的水则由于运移的差异，既可以产生板内的“水压致裂”，也会影响俯冲界面的耦合行为，进而导致慢滑移地震区的形成.由此可见，俯冲带地区深海-岩石圈流体交换及其在深部的效应是一个包含化学反应-温度-流体流动-应力变形/破坏的多物理场耦合的复杂动力学系统.俯冲带等中深部地震，虽然地震机理复杂，但多场耦合效应十分明显.需要从地球系统科学的角度出发，将流体运移、化学反应与传统的固体地球地震动力学研究相结合，着眼于多场耦合动力学综合研究，对俯冲带地区深海-岩石圈流体交换及其效应进行多时空尺度定量化表征和分析.

(3)典型区域计算网格模型、边界条件、材料特性等参数确定与快速创建.我们知道，地震研究涉及众多不同的数据.一般的断层数据是由点数据集约束而成，如何由此生成三维地质体有限元网格并保证其中断层系的几何形状、位置、交叉和整体网格质量是地震数值预测有限元模拟的先决条件，同时地震断裂模拟过程中或许还需要网格重划分来应对局部大变形和破坏所引起的网格畸变问题.因此计算网格自适应快速生成与重划分是地震数值模拟预测所面临的关键问题之一.另外，还需要综合利用地震观测、地质调查、GNSS/InSAR 形变测量、应力场观测数据及历史地震资料，对研究区诸多数据进行整合，协同并有效地使用研究区内不同方法得到的跨时空的多物理量观测数据，以快速而准确地确定有限元计算模型的边界条件、材料特性等关键参数，为后续的研究区断层系统的有限元模拟分析和未来数值预报奠定基础.

致谢感谢编辑与三位匿名评审的细致审阅、批改及建议，使得文章得到很大完善与提高．第一作者在读研期间，经导师王仲仁教授介绍认识了王仁先生，并有幸参加了王仁先生主持的材料本构关系国际会议预备会并做汇报.谨以此文纪念两位先生．