断层面的有限单元模拟方法综述①

和 平，李志雄，陆远忠，邵志刚

（1.中国地震局地震预测研究所，北京 100036；2.中国地震局地壳应力研究所，北京 100085）

0 引言

地震的孕育、发生与断层有着密切的关系［1－5］。为了更好地认识地震发生过程，前人对断层模型和机制曾进行过大量的研究。Sibson［6］在1977年提出了将断层划分为摩擦和准塑性两个域进行处理，其中摩擦域温度低，碎裂作用占主导；准塑性域温度高，主要由造岩矿物发生塑性变形所致。Strehlau［7］认为由三种构造流变层组成断层模型，即断层泥组成的上部摩擦层、复杂半脆性岩石组成的中部转换层、糜棱质剪切带组成的下部延性层。Dieterich［8］和 Gomberg等［9］提出断层滑动具有摩擦状态－速率依从机制，即断层的摩擦系数与滑动状态和历史有关。刘桂萍等［10］采用Gomberg建立的具有摩擦状态—速率依从机制、地震成核断层函数组成的数学物理模型，模拟大地震后由区域静应力场变化引起的区域触发地震现象，并探讨相关力学机制。

然而，在地球动力学、地震学、地壳形变学科的研究中，由于实际地质构造形态复杂性、地球介质分布非均匀性、以及介质物理性质多样性，基于上述模型求解偏微分方程解析解时遇到了很大的困难。20世纪60年代后期，有限单元数值模拟方法以其独特的优势——复杂几何构型适应性和各种物理问题适用性，被引入到地学研究中，并因它的优越性而得到迅速推广［11－15］。然而，在利用有限单元法在对地学问题进行研究时所面临到的第一个重要难题就是利用基于连续介质理论的有限单元法如何引入对具有非连续性间断面的处理，即对断层进行数值模拟。这一问题引起了很多学者关注，众多学者也曾对此进行过许多有益尝试，并提出过多种断层模拟方法和模型。

本文总结归纳了近年来国内外学者用有限单元法进行地学，特别是地震学、地壳形变学科研究时，处理断层所采用的主流方法。总体而言，断层处理方式可归纳为两类，一类是将断层处理为特殊连续介质体，即弱化带模型；另一类则将断层作为非连续特性的位错面来模拟，主要有劈节点模型、接触模型以及块体模型。本文就上述四大断层模拟方法及其基础上展开的地震学、形变学研究展开比较详尽的概述，以期为进行地学数值模拟时处理断层提供一些借鉴和参考。

1 介质弱化带模型

弱化带模型假定断层具有一定的厚度并被某种物质充填，断层和周围地质体的差别仅在于内在物理属性和介质参数的不同，介质彼此仍是连续的，就像没有“破”一样。弱化带模型从某种程度上说是合理的，因为根据野外地质调查，断层带内的介质常常发育成糜棱岩，断层泥等系列，和周围地质体差异较大，实验室研究也证明断层带填充物对断层的行为有重要影响［16］。

在弱化带模型中，对断层的模拟主要通过改变断层附近介质物性特征或者改变其介质参数来模拟近似断层带附近易变形的特征。部分学者对断层填充介质层用不同的本构关系去模拟，例如弹塑性模型［17－18］、应变软化模型［19－20］、粘弹性模型［21］。

王仁等［17－18］用二维有限单元法研究华北地区地震迁移规律时，对断裂带的处理是将相邻的同方向的断层合并在一起作为介质弱化带，宽度10km左右。考虑到断裂带内岩石体是经过错动的，因此将弹性模量E和摩擦系数μ取低；又由于断裂带取得较宽，E，μ的降低量不会很大。考虑的破裂有两种，即沿断层方向的剪破裂和拉破裂。其中采用的剪破裂准则为库仑破裂准则，当断层发生错动时，断层内应力这时与应变之间服从理想塑性体法则。模拟结果基本重现出了近十二年来的地震迁移规律，得到了华北地区应力场，并为未来地震危险区提出了参考。

殷有泉等［19］从岩石弹塑性一般本构方程出发，考虑断层介质摩擦系数和内聚力因变形而降低的特点，建立了断层介质的非稳定本构关系，并用一个能量表述的非稳定准则来判断包含断层和围岩在内的整个系统的非稳定性的到来。在此基础上模拟出了1976年唐山7.8级地震弹性回跳和1971年圣费尔南多6.6级地震，得到的模拟结果与当时观测的地震位移数据在规律和数量级上大体一致。

王继存等［20］在应变软化模型基础上研究了唐山地震断层破坏及其构造应力场，探讨了唐山地震从孕育到发震的力学过程，分析了唐山发震断层破坏的力学机制，并利用增量分析的方法，按时间序列近似地模拟出了发震断层和周围岩石从孕震到发震的过程。模拟中，唐山地区断层带首先进入塑性阶段，并在闭锁段出现应力集中，相当于唐山断层处于临震状态；然后在闭锁段两侧施以方向相反的指定位移，模拟发震时断层的突发走滑，相当于断层闭锁段被剪破裂连通而发震。通过调整介质内摩擦系数和内聚力的方法来模拟震源区的应力释放，模拟得到的塑性区分布与各时段震余震分布有着较好的吻合。

Mian Liu等［21］用三维粘弹性模型研究了区域压应力场背景下青藏高原扩张运动的机理，并通过调整欧亚板块与印度板块间的断裂带介质粘弹性来近似模拟出两板块间的板块挤压作用、重力侧向力作用机制以及地面速度变化。研究结果表明，若青藏高原高度只为本身高程百分之五十时，则该区域走滑断层与逆断层为主要断层，且东西向无扩张运动趋势。当今高原地壳扩张状态只有当青藏高原达当今75%高程时才会出现。

陆远忠等［22］建立了含有母体岩石、硬包体和随机分布的小裂纹的三维有限元模型，计算了包体和各层实体中的应力分布，模拟强震前岩石的破裂和小震的空间分布特征。结果表明，强震前在孕震体即包体附近均出现了高应力集中单元，它们是形成小震空区、条带和地震空间丛集图像的基础。随机裂纹的存在有利于应力在孕震体（包体）外的裂纹端部集中，发生小震，形成包围孕震包体的前兆地震活动图像。而包体中的应力逐渐增加，为发生强震提供了条件。

陈连旺等［23－24］在研究华北地区构造应力场及其演化时，根据华北地区地壳波速三维结构构建了纵向分层、横向分区的三维有限元模型，不同层、分区内的物性参数不同，对断层带的处理主要是增加其介质泊松比和降低杨氏模量来加以实现。计算结果表明，华北地区大多数区域，最大主压应力与最小主压应力（最大主张应力）均位于接近水平的平面内；中等主应力位于接近垂直的平面内即处于走滑应力状态。之后，陈连旺等［25］又在上述弱化带处理的基础上建立了三维粘弹性有限元模型，同时引入了区域背景应力场和库仑破裂准则，模拟了1966年邢台地震引起的华北地区应力场动态演化过程，得到了强地震对华北地区应力场的震时扰动作用及其引起的一年时间尺度的库仑破裂应力的动态变化速率和对下一次强震触发作用。

陈化然等［26］建立了包括上地壳、下地壳、上地幔等4层结构的三维有限元模型，计算了川滇地区背景应力场、断层蠕动和强震触发的应力场及它们的动态变化。结果表明：川滇地区后续地震大多发生在前面地震引发的库仑破裂应力正值区，前一地震对后续地震有一定的触发作用，强震是在较高的应力背景下成组发生的。

李平等［27］用三维有限元数值模拟研究了辽宁及邻近地壳构造应力场及其与地震活动的关系，模型考虑了地壳上地幔三维地震波速结构、地质构造及地壳深部构造环境等，其中断层以宽约1km的条带模拟，并结合震源机制解与实际地应力测量资料确定了随深度变化的水平边界力。研究表明，地质构造、断层分布与介质力学参数直接影响到应力场分布，因此对这些参数研究的精细程度将直接影响到应力场分布的可靠程度。

总体而言，弱化带模型能很好地描述出了断层带内介质的属性特征以及填充物对断层行为的影响，且模拟出断层带更容易发生变形的特点。然而其不足之处在于模型计算过程中断层上单元的共同节点有且仅有一个方向的位移，因此该模型无法模拟出断层和周围介质的交界发生的位错，更难模拟出断层上存在的摩擦机制。

2 劈节点模型

劈节点（split node technique）模型由 Melosh等［28－29］在1981年正式提出，是对Jungels等［30－31］拟断层滑动方法的一种修正和改进。其核心思想是通过引入一个新载荷向量，使得模型中断层处所对应的节点在毗邻的单元上产生大小相同而方向相反的位移即位错，就像断层面上的节点被“劈”开了一样，断层则是由这些“劈”节点所组成。

劈节点模型主要是通过初始位移的形式将断层位错引入有限单元数值模拟中。以一维的两个单元为例，以下是该模型的处理思路。单元1，2分别位于断层的左右两侧，U为位移，上标表示单元号，下标表示单元内的节点号；ΔU表示位错量；则有

其中分裂节点上位移平均值为¯U12＝¯U21，ΔU12＝－ΔU21。其示意如图1。

图1 一维劈节点模型示意图Fig.1 Sketch of One Dimensional Split Node Model.

利用以上关系，可以将单元局部刚度矩阵合成以下全局刚度矩阵

其中Ui为位移；Fi为位移；i为全局节点号。

Yoshioka等［32］通过横向异性的三维有限元模型模拟了1993年Kushiro－oki地震的同震和震后位移以及应力场的变化，模拟中断层位错量大小通过近场地震波反演来给出，研究认为同震形变在断层的上下面较大而震后形变主要位于俯冲面与周围软流圈附近，且同震状态下应力将在未来20年内保持不变。

Hearn等［33］研究了震后余滑作用和粘弹性松弛所导致的震后形变的区别，以及如何利用这些区别来辨别出地壳内部及地幔的分布和流变特性。研究认为一个设计良好的GPS流动测点观测台网可足够将震后余滑同线性粘弹性松弛区分开来，同时还可得到粘弹性层的厚度和粘性参数。研究指出为最大限度利用GPS观测资料，部分GPS站点应沿地表破裂带布设。

Sheu等［34］建立了三维有限元模型对1999年台湾集集7.6级地震震后形变的机理进行了研究，通过对比分析集集地震后97天内的GPS资料，认为集集地震震后形变是震后余滑和应力松弛共同引起的，其中余滑起主要作用。Buiter等［35］等用二维数值模型研究了在百万年时间尺度下，不同俯冲带板系统（重力拖拽力、板块运动速度、俯冲带交界面摩擦系数等）对板块接触边缘带垂直位移的影响。结果表明板块接触边缘带垂直位移主要受俯冲断层摩擦系数和倾角参数等影响，对板块边缘垂直影响幅值介于2km到4km之间。

Hu等［36］用Maxwell粘弹性体模型研究了1960年智利地震的震后形变，探讨了介质的流变性对地壳形变的影响，其中俯冲板块，海洋及大陆板块为弹性介质，大陆及海洋地幔为粘弹性介质。模型中以抛物线函数拟合了俯冲板块剖面，由二维延展至三维，以赋予断层面一个位错量来模拟地震的发生。模拟结果认为1960年智利地震35年后GPS观测得到地面水平位移特征可能是由地幔应力松弛导致的。

Hyodo等［37］进行了日本中部 Niigata－Kobe构造带上震间应变积累研究。研究表明当板间地震周期比较长时，Niigata－Kobe构造地壳相对较薄的部位存在着较高的应变积累。Zhao等［38］用三维模型探讨了理论上三类断层中地壳分层结构对地表形变和应力的影响，结果表明地壳分层与各项同性地壳结构可导致至少有20%的地表形变差异和15%的应力差异，因此当在进行数值模拟时必须考虑地壳纵向分层、横向分块的结构对模拟结果的影响。

Dixon等［39］研究了加州Baja半岛的断层滑动率受地震轮回和岩石圈流变性的影响，并对当地新构造活动和地震危险性做出了评估。模拟结果表明Agua Blanca和San Miguel断层均在一定程度上受上述两个因素的影响均比较明显，并导致约2 mm／yr断层活动率的变化。此外，Dixon通过计算认为Agua Blanca断层目前可能处于地震轮回的晚期阶段，同时具备着每隔百年孕育6.1至7.0级地震的能力。

Aagaard等［40］利用 Andrews等［41］提出的拉拽劈节点（TSN）断层处理方法研究了在同质和分层半空间的情况下走滑和逆冲断层的动态断裂过程和由此引起的地震动，研究表明介质摩擦系数，剪切模量以及断层的深度对破裂过程有着一定的影响，且与断裂的速度和传播方向有着直接的联系。

劈节点模型简单有效地以初始位移的形式将断层近似引入有限单元数值模拟中，同时通过引入一个载荷向量模拟地震引起的位错。然而，由于这种模型只是通过模拟断层的位错行为来描述断层的存在，从本质上来说并没有对客观原本存在的断层进行模拟。因此劈节点模型在地学中的应用主要局限于地震同震和震后效应等相关研究。

3 接触单元模型

接触模型主要是通过引入一对接触面，并在这对接触面的表面“焊接”上特殊面单元（零厚度的接触单元）来模拟地质体中断层的存在和其一定摩擦机制下的滑动特征的。实质是将断层视为不同地质体之间所接触的部分，断层的活动则通过不同地质体间的接触滑动来模拟。

接触模型的分析过程比较复杂，有限元方程组无法通过一次求解而实现，而须采取增量求解法对于每个增量必须使用“试探—校核”的迭代方法进行求解。首先根据前一步结果和本步给定的载荷条件，并通过接触条件的检查和搜寻，设定此步第一次迭代求解时的接触面的区域和状态，然后再根据上述接触面区域及状态的假设对于接触面上的每一点将接触界面约束条件改为等式约束作为定解条件引入方程组并进行方程的求解，最后利用接触面上于等式约束对应的不等式约束条件作为校核条件对解的结果进行检查。如果两者相符则转入下一增量步计算，如两者不相符合则重新进行搜寻及迭代求解，直至满足校核条件而转入下一增量步计算。

通过零厚度的单元来模拟非连续面的行为的思想最早由Goodman等1968年［42］提出，起初为节理单元法，后经过 Hallquist、Cescotto、Peric等多人［43－45］陆续调整修改逐渐形成了一套较完善的接触—碰撞理论，并被广泛应用于岩土工程、建筑等各个领域，后在地学中也得到了越来越多的应用。

Tom Parsons［46］最早利用含有接触面单元的三维粘弹性模型模拟了1906年旧金山7.8级地震后圣安地列斯断裂带的应力恢复过程，模型中的间断面代表旧金山断裂带上三条的走滑断裂。断层面通过零厚度接触面单元来模拟，接触单元焊接在周围粘弹性单元的表面上，断层可发生形变且遵循库伦摩擦准则。此外Parsons等［47］还利用含地表高程，不同厚度地壳层和蠕滑断层的复杂三维模型并结合GPS观测资料研究了加利福尼亚地区的构造应力情况，研究发现部分地壳运动并非板块构造所驱动，而是由于地形下降所导致，Garlock地震并非由于高构造应力导致，而是与圣安地列斯断层的滑动引起的。

Masterlark［48］以1985年墨西哥Jalisco－Colima的Mw8.0俯冲带地震为例，结合GPS观测资料，通过对比位错模型解析解与三维有限元数值计算结果，研究了位错模型中均匀各项同性半空间泊松体（HIPSHS）中的每一种假设对地形变的影响。研究表明，各假设对地形变模拟影响敏感度由低到高分别为半空间、泊松体、弹性各向同性和同质体，且除半空间外，各假设所导致的敏感度大小远大于GPS数据的不确定度。

Ellis等［49］用含有摩擦机制的断层面模拟了新西兰阿尔卑斯断裂带的弹性和非弹性应变积累，地壳体中断层滑动机制和介质行为由岩石实验和地球物理测量所确定。研究认为发生地震的粘滑断层之下存在着一个蠕动延展带，且控制着地面水平位移，所模拟得到的地面水平速度场与GPS观测结果取得了很好的一致性。

Ganas等［50］模拟了万年尺度下希腊 Hellenic Arc地区地形变、强地震活动性以及Crete地区隆起的机理，模型表面由GPS地面水平位移约束。研究表明Aegean地形变、Crete隆起主要由一个33 mm／年快速移动的上冲板块与一个相对静止（大概5mm／年）下俯冲板块导致的，因此Hellenic Arc的构造更接近于是陆地俯冲板块。

刘峡等［51］在研究华北地区现今地壳运动及形变动力学数值模拟时分别采用了三维连续弱化带模型和二维滑动摩擦模型。在二维模型中采用了接触单元法来实现对非连续变形断层的模拟，其相对运动由接触面的正压力、摩擦系数和抗摩擦强度来控制，一旦接触面的剪切力超过摩擦强度或静摩擦力，则相对运动发生。研究表明，接触模型模拟结果与观测到的GPS数据平均离散度高于三维介质弱化模型，这可能是由于接触模型忽略了华北地区地壳纵向不均匀性导致的。

李红等［52］根据广州地区地震地质、地球物理、震源机制解以及GPS观测等相关资料，用三维有限元模型模拟计算了广州地区四条断裂的现今活动方式以及库仑破裂应力年累积速率。其中，对广州地区的主要活动断裂主要采用接触单元表示的非连续变形的断层间断面来模拟，间断面之间遵从库仑滑动摩擦准则，模型单元的杨氏模量和泊松比依据广州地区波速结构确定。模拟结果表明，广州—从化断裂和瘦狗岭断裂东段断裂面上的库仑破裂应力年累积速率明显高于瘦狗岭断裂西段和珠江口断裂西支、东支，数值模拟计算结果与该地区地震活动特征基本一致。

王凯英［53］对川滇地区应力场形成机制及断层相互作用进行数值模拟，模型中考虑了主要活动断裂的存在，并将断层处理为接触摩擦面，且断层带由具有一定宽度的相对软弱介质组成。分析表明川滇地区断层相互作用现象是块体非均匀运动过程中应力场调整的反映，是块体运动的结果。此外王凯英［54］以青藏高原北部地区几次强震为例研究了介质非连续的断—块构造模型中地震永久位移造成的区域断层静库仑应力变化。研究表明在断块构造模型中走滑型为主的地震位错引起的区域断层库仑应力变化在地震破裂面以外随距离衰减很快，影响范围仅为几十公里，和余震分布尺度相当；对青藏高原北部的多次地震的模拟结果显示，先前地震的断层位移在后续地震断层面所造成的库仑应力甚微，预示了远距离地震之间的静应力触发效应不明显。

接触模型通过引入一种特殊的零厚度单元来进行模拟断层存在，断层被近似为由接触单元对所组成的接触面。这种模型最大的优点在于可以充分考虑断层面的摩擦机制对模拟结果的影响。不足之处在于接触模型的计算分析过程十分复杂，而且接触刚度的选取困难。为避免接触面互相穿透而影响计算精度，应选取较大的接触刚度，然而太大的接触刚度会产生计算收敛困难，接触表面互相跳开等。此外接触模型也没有考虑断层带介质参数及断层内填充物（断层泥）对模拟结果的影响。

4 块体模型

块体模型是指将地学问题看成一个多地块相互作用的动力学过程，各块体本身同时具有滑移和变形的性质，模型中的断层实质就是块体系统中的地块边界。通过多块体作用体系来近似模拟地学现象可能最早起源于活动亚板块和构造块体相对运动假说，而这一思想能被成功应用于有限元数值模拟更大程度上归功于美籍华人石根华提出的非连续变形分 析 方 法 （Discontinuous Deformation Analysis）［55－56］。

石根华提出的DDA方法假定每个块体处处具有常应力，常应变，块体的运动和变形由12个独立的变形参数确定，如下：

［D］＝ ［u0，v0，w0，α0，β0，γ0，εx，εy，εz，γxy，γyz，γzx］T（3）其中 （u0，v0，w0）指块体质心（x0，y0，z0）的刚体平移；（α0，β0，γ0）指块体绕x，y，z轴的转角；（εx，εy，εz，γxy，γyz，γzx）指块体的3个正应变和3个剪应变。

由于式（3）的常应变假定限制了块体复杂变形的能力，为更好地求解块体内部的位移场和应力分布，DDA的块体单元被离散成mb个四面体单元，内嵌有有限元网格的块体单元势能为

式中，［ke］为单元的刚度矩阵；［pe］为单元等效节点载荷矩阵；［ae］为节点位移。

对于离散模型，块体单元间的接触界面是位移间断面，界面的滑动和分离是由Mohr－Coulomb准则和无拉应力准则控制的。DDA采用罚函数求解接触问题，等价于在接触点施加刚度很大的法向和切向弹簧，弹簧的加设和移开称之为开—闭。迭代假设所定义的块体系统有n个块体和k个接触，则块体系统的总势能包括各块体单元势能及接触弹簧的应变能：

式中，［kcij］为接触弹簧矩阵，根据变分原理δ∏p＝0，经推导可得块体系统的总体平衡方程：［K］［a］＝［P］，其中［K］为块体系统整体刚度矩阵，由有限单元刚度矩阵和接触矩阵集合而成；［P］为块体系统整体荷载列阵，由有限元等效节点荷载列阵集合而成。鉴于块体模型以块体为研究重要组成部分，因此块体模型主要研究大尺度的的地学问题。

Liu［57］基于华北地区的地质构造的断块特征，用非连续变形方法模拟了该地区的长期形变，包括无震断层滑动和断层形变。研究发现对一条断层，由相邻断块内块体旋转引起的断层滑动方向可能与区域构造挤压引起的滑动方向相反，但由于由构造挤压引起的滑动量级远大于由断块旋转引起的断层滑动，因此，断层滑动方向基本与一个地区构造挤压方向一致。

蔡永恩等［58］用新的 LDDA（Lagrangian Discontinuous Deformation Analysis）方法模拟了唐山断层破裂、错动和应力释放的整个过程。LDDA采用了DDA方法中的接触判断准则，采用区域分解方法求解出接触面上满足的约束条件边界力，再以接触力为边界条件求出每个块体内部的位移和应力。研究发现唐山地震断层的破裂速度和应力将与断层上的初始剪应力有关，且唐山发震断层的最大动态、准静态位错和剪应力降均发生在中间部位，断层破裂的传播速度从震中向东南和西北方向分别为3.08km／s和1.18km／s。

白武明等［59］用DDA＋FEM模型研究了华北地区各地块相互制约的构造环境中发生1976年唐山地震的动力学过程，并探讨了大震引起的华北地区各地块特别是鄂尔多斯地块运动变形及各地块边界断层上应力状态变化的特征。研究认为唐山地震的发生引起了华北地区各地块边界断层上的应力状态发生一定程度的变化，其中通过唐山震区的张家口—蓬莱北西走向的断裂带及鄂尔多斯地块北部边界断层的主要地带上剪应力增加而法向应力减小，使得这些断裂带的地震危险性增加，这与唐山地震后鄂尔多斯地块北部边界断层发生一系列6级强震事实基本相符。

张瑞青等［60］通过对1999年岫岩5.4级地震前震及余震分布图像研究和前人对海城地震的研究，提出了海城、岫岩地震发震构造模型，并利用DDA＋FEM方法模拟了2次地震释放的主应力变化、最大剪应力变化等值线图，地震前后位移变化矢量图，及发震断层滑移随时间的变化。模拟结果分别与相应地震的震源机制、宏观等震线、发震断层的走向性质等基本一致。

陈祖安等［61］用三维DDA＋FEM 方法以GPS资料做位移速率边界和震源机制约束，模拟了1997年玛尼7.9级大震发生对块体边界断层应力状态的影响。此外，陈祖安［62］还研究了2008年汶川8.0级地震孕震机理。研究发现，龙门山断裂带东西两侧地势、地壳厚度、分层以及物性明显变化对汶川大地震的孕育和发生均起了关键的作用。其中高原地壳物质向东水平运动，受到龙门山断裂带东侧介质刚性较大的四川盆地阻挡，使得汶川大地震发震断层在大震前已积累了相当水平的应变能，处于不稳定状态。

块体模型将刚体位移和块体变形采用统一的有限元格式求解，所有块体不仅允许自身有位移和变形，而且允许块体间有滑动、转动、张开等运动形式，能够得到块体系统的大位移和大变形解。块体模型结合有限元法和离散元法的优点，能够较好的反映出块体的运动趋势，因此在地学中取得了广泛应用。然而，在块体模型中对断层的模拟实质是作为次级块体的边界带，故在二维和三维模型中大多数断层往往被简化为直立断层，同时也无法考虑块体内部断层存在和分布对模拟结果的影响。

5 小结与展望

本文概述了近30年来国内外学者用有限单元法研究地学问题时就具有非连续性的断层所采用的几大主流处理方法，并扼要叙述了在各个断层模型基础上展开的地震学和断层动力学方面的研究。

总体来讲，有限单元法关于断层的处理主要有两大分流：①20世纪80年代以王仁等人为代表提出的连续介质弱化带模型，主要通过改变断层带物性特征或其介质参数来得到断层带易变形的特征。这种方法在一定程度上是科学的，因为它考虑到了断层带介质与周围岩体性质的差异性，改变本构方程或介质参数，使其在一定程度上更容易发生变形来近似模拟断层带附近易变形的特征；②20世纪60年代以来针对处理断层非连续特点所提出的非连续介质间断面模型。主要代表有1968年Goodman等人提出节理单元、1981年Melosh提出的劈节点处理方法、1988年石根华提出的非连续变形分析方法、接触单元模型。这些模型的共同点均是将断层处理为非连续间断面，其与弱化带法相比最直接的优点在于能够更好的模拟出断层两侧区域存在的速度和位移场的跃变性以及断层受到摩擦滑动机制，可相对更真实的表现客观断层活动特点，然而，却很难体现出断层带内介质断层泥等对断层滑动的影响特征等。

以上断层处理方法在地学问题模拟分析中得到了广泛的应用，在对断层进行处理时，每一类模拟方法也均存在着它自身的优势和局限性。为更准确地认识断层与地壳形变及地表应力分布关系、深化断层活动与地震孕育的关系研究，并为未来地震预测提供更加真实可靠的依据，本文提出今后对断层的模拟方法可能将会有以下几个发展趋势：

（1）在同一研究区域中，断层处理手段的多元化。即根据不同地质构造选取不同的处理模型，如新发生的断层面易用接触单元来模拟，有断层泥发育的断裂可采用介质弱化模型；根据不同的研究区域尺度选取不同的模型，如研究区域为中国大陆时宜先采用块体模型进行处理，再对块体内部的断层采用接触或弱化带模型。

（2）同一断层面上分层化处理。如将断层面由浅至深划分为摩擦和准塑性两个不同的域，其中摩擦域用接触模型来模拟，准塑性域则采用应变软化法来使得岩矿物能够塑性变形；或可将同一断层处理为三个构造和流变层组成的断层模型，上部摩擦层、中部转换层和下部延性层，以探讨各断层模型对构造应力场和形变场的影响。

（3）过去研究大尺度区域的地学问题时，相对小的断裂带往往都被视为单一平面。然而野外实际调查表明即使是规模相对较小的断裂带其产状分布也往往不是一简单平面，断层的不同部分倾角、形态、深度均可能会存在较大差别。已有研究表明，断层产状和分布的选取在地学中起着重要的角色。因此随着地球物理勘探技术、地质力学、地球化学和观测技术的发展，在对地学问题进行数值模拟时，断层产状、分布及其介质参数的选取也必将会趋近真实化。

（4）断层带物性实测和模拟方法研究。根据不同条件下的大量实验结果［16］，断层泥可大概分为三类：碎屑类、碳酸盐类、层状硅酸盐类。断层泥的矿物成分又很大程度地影响着断层的力学行为。因此在对断层带进行数值模拟时，应结合考虑实际的场地构造环境因素如温度、围压、孔隙压以及断层泥矿物成分类型来对断层的摩擦滑动特征、滑动方式和断层泥变形机制等加以模拟。

［1］邓起东.活动断裂研究的进展与方向［M］.北京：地震出版社，1991.

［2］闻学泽.活动断裂地震潜势的定量评价［M］.北京：地震出版社，1996.

［3］薄万举.地壳形变与地震预测［M］.北京：地震出版社，2001.

［4］李兴才.断层的震前滑动对唐山地震发生的影响［J］.地震学报，1992，14（3）：304－308.

［5］周硕愚.断层形变测量与地震预报［J］.地壳形变与地震，1994，14（4）：90－97.

［6］Sibson R H.Fault rocks and fault mechanisms［J］.J.Geol.，1977，133：191－213.

［7］Strehlau J.A discussion of the depth extent of rupture in large continental earthquakes［J］.Geophys.Monogr.，1986，37：131－145.

［8］Dieterich J H.Earthquake on faults with rate and state－de－pendent friction［J］.Tectonophysics，1992，11：115－134.

［9］Gomberg J.Beeler N M.Earthquake triggering by transient and static deformation［J］.J.Geophys.Res.，1998，103（B10）：24411－24426.

［10］刘桂萍，傅征祥.摩擦状态－速率依从的区域地震触发模型研究［J］.地震，2004，24（1）：176－183.

［11］王仁.有限单元等数值方法在我国地球科学中的应用和发展［J］.地球物理学报，1994，37（增刊）：128－138.

［12］王妙月，底青云，张美根.地震孕育、发生、发展动态过程的三维有限元数值模拟［J］.地球物理学报，1999，42（2）：218－227.

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

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

［15］张彬，杨选辉，陆远忠.地震动态应力触发研究进展［J］.西北地震学报，2008，30（3）：298－303.

［16］马胜利.模拟断层带摩擦滑动性状与变形特征［J］.中国地震，1986，2（2）：72－78.

［17］王仁，何国琦，殷有泉.华北地区地震迁移规律的数学模拟［J］.地震学报，1980，2（1）：32－42.

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

［19］殷有泉，张宏.模拟地震的应变软化的数学模型地球物理学报［J］.1982，25（5）：414－423.

［20］王继存，续春荣.唐山地震断层破坏及其力学过程的数值模拟地震地质［J］.1989，11（4）：71－76.

［21］Liu M，Yang Y.Extensional collapse of the Tibetan Plateau：Results of three－dimensional finite element modeling［J］.J.Geophys.Res.，2003，108 （8）：2361.DOI：10.1029／2002JB002248.

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

［23］陈连旺，陆远忠，张杰，等.华北地区三维构造应力场［J］.地震学报，1999，21（2）：140－149.

［24］陈连旺，陆远忠，郭若眉，等.华北地区断层运动与三维构造应力场的演化［J］.地震学报，2001，23（4）：249－361.

［25］陈连旺，陆远忠，刘杰，等.1966年邢台地震引起的华北地区应力场动态演化过程的三维粘弹性模拟［J］.地震学报，2001，23（5）：480－491.

［26］陈化然，陈连旺，马宏生，等.川滇地区应力场演化与强震间相互作用的三维有限元模拟［J］.地震学报，2004，26（6）：567－575.

［27］李平，卢良玉，卢造勋，等.辽宁及邻区地壳构造应力场及其与地震活动关系的三维有限元数值模拟研究［J］，地震学报，2001，23（1）：24－35.

［28］Melosh H J，Raefsky A.A simple and efficient method for introducing faults into finite element computations［J］.Bull.Seism.Soc.Am.，1981，71（5）：1391－1400.

［29］Melosh H J，Williams C.A Mechanics of graben formation in crustal rocks：A finite element analysis［J］.J.Geophys.Res.，1989，94（B10）：13961－13973.

［30］Jungels P H.Models of tectonic processes associated with earthquakes［D］.Pasadena：California Institute Technology，California，1973.

［31］Jungels P H，Frazier G A.Finite element analysis of the residual displacements for an earthquake rupture：source parameters for the San Fernando earthquake［J］.J.Geophys.Res.1973，78：5062－5083.

［32］Yoshioka S，Tokunaga Y O.Numerical Simulation of Displacement and Stress Fields Associated with the 1993Kushiro－oki，Japan，Earthquake［J］.Pure appl.geophys.，1998，152：443－464.

［33］Hearn E H.What can GPS data tell us about the dynamics of post－seismic deformation［J］.Geophys.J.Int.，2003，155：753－777.

［34］Sheu S Y，Shieh C F.Viscoelastic－afterslip concurrence：a possible mechanism in the early post－seismic deformation of the Mw7.6，1999Chi－Chi（Taiwan）earthquake［J］.Geophys.J.Int.，2004，159（3）：1112－1124.

［35］Buiter J H，Govers Rob，Wortel M J.A modeling study of vertical displacements at convergent plate margins［J］.Geophys.J.Int.2001（147）：415－427.

［36］Hu Y，et al.Three－dimensionalviscoelastic finite element model for postseismic deformation of the great 1960Chilean earthquake［J］.J.Geophys.Res.，2004，109.DOI：10.1029／2004JB003163.

［37］Hyodo M K，Hirahara.A viscoelastic model of interseismic strain concentration in Niigata－Kobe Tectonic Zone of central Japan［J］.Earth Planets Space，2003，55：667－675.

［38］Zhao S，et al.3－D finite element modelling of deformation and stress associated with faulting：effect of inhomogeneous crustal structures［J］.Geophys.J.Int.，2004，157：629－644.

［39］Dixon T J，Decaix F，Farina K.Seismic cycle and rheological effects on estimation of present－day slip rates for the Agua Blanca and San Miguel－Vallecitos faults，northern Baja California，Mexico［J］.J.Geophys.Res.，2002，107（B10），2226，DOI：10.1029／2000JB000099.

［40］Aagaard B T，Hall J F，Heaton T H.Characterization of near－source ground motions with earthquake simulations［J］.Earthquake Spectra，2001，17（2）：177－207.

［41］Andrews D J.Test of two methods for faulting in finite－difference calculations［J］.Bull.Seism.Soc.Am.，1999，89：931－937.

［42］Goodman R E，Taylor R L，Brekke T L.A model for the mechanics of jointed rock［J］.Soil Mech.and Found Proc，1968，94：637－658.

［43］Hallquist D，John O.LS－DYNA3DTheoretical Manual［M］.Livermore：Software Technology Corporation，1994.

［44］Cescotto S，Charilier R.Frictional Contact Finite Elements Based on Mixed Variational Principles［J］.International Journal for Numercial Method in Engineering，1992，（36）：1681－1701.

［45］Peric D，Owen D R.Computational Model for 3－D Contact Problems with Friction Based on the Penalty Method［J］.International Journal for Numercial Method in Engineering，1992，35：1289－1309.

［46］Parsons T.Post－1906stress recovery of the San Andreas fault system calculated from three－dimensional finite element analysis［J］.J.Geophys.Res.，2002，107（B8）：2162，DOI：10.1029／2001JB001051.

［47］Parsons T.Tectonic stressing in California modeled from GPS observations［J］.J.Geophys.Res.，2006，111（B03）407，DOI：10.1029／2005JB003946.

［48］Masterlark T.Finite element model predictions of static deformation from dislocation sources in a subduction zone：Sensitivities to homogeneous，isotropic，Poisson－solid，and halfspace assumptions［J］.J.Geophys.Res.，2003，108（B11）：2540，DOI：10.1029／2002JB002296，.

［49］Ellis S，et al.Simplified models of theAlpine Fault seismic cycle：Stress transfer in the mid－crust［J］.Geophy.J.Int.，2006，DOI：10.1111／j.1365－246X.2006.02917.x.

［50］Ganas A，Parsons T.Three－dimensional model of Hellenic Arc deformation and origin of the Cretan uplift［J］.J.Geophys.Res.，2009，114 （B06）：404，DOI：10.1029／2008JB005599.

［51］刘峡，傅容珊，杨国华.用GPS资料研究华北地区形变场和构造应力场［J］.大地测量与地球动力学，2006，26（3）：33－39.

［52］李红，陈连旺，李红江.广州地区活动断裂的数值模拟［J］.大地测量与地球动力学，2008，28（2）：39－44.

［53］王凯英，马瑾.川滇地区断层相互作用的地震活动证据及有限元模拟［J］.地震地质，2004，26（2）：259－272.

［54］王凯英.断块模型中走滑型地震应力触发研究——以青藏高原北部几次强震为例［J］.地球物理学报，2009，52（7）：1776－1781.

［55］Shi G H.Discontinuous deform ation analysis：a new numerical model for the statics and dynamics of block systems［D］.Berkeley：Department of Cilvil Engineering，University of California，1988.

［56］白武明，林邦慧，陈祖安.1976年唐山大震发生对华北地区各地块运动与变形影响的数值模拟研究［J］.中国科学（D辑），2003，33（增刊）：99－107.

［57］Liu L.Modeling aseismic fault slips and block deformation in northern china by DDA［A］∥Proc of the First International Forumon Discontinuous Deformation Analysis（DDA）and Simulations of Discontinuous Media［C］.Albuquerque：TSI Press，1996：373－383.

［58］自武明，林邦慧，陈祖安.1976年唐山大震发生对华北地区各地块运动与变形影响的数值模拟研究［J］.中国科学（D辑），2003，33（增刊）：99－107.

［59］蔡永恩，何涛，王 仁.1976年唐山地震震源动力过程的数值模拟［J］.地震学报，1999，21（5）：469－477.

［60］张瑞青，魏富胜，乔成斌.用（DDA＋FEM）方法数值模拟1975年海城、1999年岫岩地震发生的过程［J］.地震学报，2005，27（3）：163－170.

［61］陈祖安，林邦慧，白武明，等.1997年玛尼地震对青藏川滇地区构造块体系统稳定性影响的三维DDA＋FEM方法数值模拟［J］.地球物理学报，2008，51（5）：1422－1430.

［62］陈祖安，林邦慧，白武明，等.2008年汶川8.0级地震孕震机理研究［J］.地球物理学报，2009，52（2）：408－417.