2011年3月11日，日本发生M9.0地震，引发海啸、核泄漏等一系列次生灾害，造成大量人员伤亡.已有不少学者研究其对周边地区的影响(Asano et al., 2011; Enescu et al., 2012; Hiratsuka and Sato， 2011; Ozawa et al., 2011; Takahashi，2011). 如此规模的地震对远场的影响同样不可忽略(Gonzalez-Huizar et al., 2012; Sun and Zhou， 2012; Wang et al., 2011; Yukutake et al., 2011; Zhao et al., 2012; 陈为涛等，2012)，本文将计算并讨论此次地震之后，我国华北地区应力状态有何变化，以及华北地区各大断层稳定性受到何种影响.

地球作为一个复杂的物理系统，在受到力的作用时，短期内表现为弹性体.因此计算地震发生的同震效应，宜将地球作为弹性体处理，使用弹性位错理论.20世纪50年代以来一些学者发展了这套理论，Okada(1992)给出了半无限空间位错解析公式，汪荣江(2003)给出了分层介质中位错的计算程序，孙文科等(2009)给出计算球状分层地球的解析公式.虽然解析解有计算快捷简便的优点，但是对于存在地形和Moho面起伏、横向不均匀性等复杂情况的模型却无法计算，使得该方法的通用性受到限制，数值计算的方法可以避免这方面的困难.

有限元法正是求解偏微分方程的一种重要的数值方法.但是传统的有限元法使用连续的形函数，因 此不能处理这种有间断解的问题.为了克服这类困难，人们提出了一些特殊的处理方法，如DG(Discontinuous Galerkin，不连续伽辽金)方法、罚函数方法、双节点方法等(Arnold，1982; Oden et al., 1982; 朱桂芝和王庆良，2005).这类方法的特点是，单元内部的形函数仍然连续，但是单元之间的结合处可以出现间断，因此可以处理间断解的问题，额外的代价就是单元之间的分界面必须很好地贴合间断面，增加了建模的困难，尤其当间断面形状不规则时，会造成网格质量较差，影响计算精度.Belytschko等人提出了一种扩展有限元方法(Belytschko and Black， 1999; Belytschko et al., 2001; Dolbow et al., 2000; Gracie et al., 2007; Sukumar et al., 2000)，允许单元内部出现间断面，从而避免建模的困难.其方法是仍然使用连续的形函数，但在节点自由度之外引入额外的未知量来刻画间断面，从而允许单元内部出现间断，额 外的代价是增加了未知量的数目，单元内部出现间断对数值积分也造成了麻烦(Belytschko et al., 2001).

本文通过另一种途径来计算弹性位错问题.将位错等效为体力，添加到平衡方程的右端项中.只需 要计算出这个体力的大小，便将位错问题转化为传统的有限元问题，完全避免上述各种方法为了处理 不连续性而遇到的困难.在实际使用中，计算这个体力项所需的代码量小，添加体力项也是各种主流有限元计算程序都支持的功能，对计算网格更是没有任何额外的要求，并可以处理任意复杂形状的断层，十分方便.

为检验此方法以及程序的有效性，通过一个简单算例并将结果与解析解进行对比；最后计算并讨论2011日本东北M9.0地震对华北地区主要断层稳定性的影响.

本节公式均使用了Einstein求和约定.

Burridge和Knopoff(1964)指出，弹性介质中的位错源可以与体力源完全等价，并给出关系式：

该式将任意位错等效为计算区域中各r′点处的体力，右端积分号中的第一项表示与位错等效的体力，第二项表示与位错的一阶导数即应力错等效的体力.

定义两个张量 A 、 B 的双点乘规则如下：

则(1)式可写为矢量式：

其中

假设待求解的有限元问题的弱形式为： a(u,v)=(f,v)，将上述等效体力写到方程右端体力项 f 之中.然后在单元积分时，将上述 f 的表达式与试探函数 v 做内积，为此只需将 f 与各单元形函数做内积：

其中 Φ _m表示单元K中的第m个形函数，V^K表示单元K的体积微元.将(4)式右端第一项变形并应用Green公式：

其中S^K表示单元K的外表面， n(r)表示r点处单元K外表面的法方向，若r不在单元表面，则取 n(r)=0.如果S^K∩Σ为空集，即单元表面与断层 面没有公共部分，则对断层面上任意一点r，有 n(r)=0，(5)式第一项积分结果为0；如果S^K∩Σ不为空集，单元表面与断层面有公共部分，注意到最后需要对所有单元积分结果进行求和，相邻单元在交界面上的法方向相反，(5)式第一项积分结果互相抵消.因此无论断层面与单元表面相交与否，(5)式第一项积分的贡献始终是0，只需计算第二项积分.

同理，(4)式的第二项可化为：

本文的计算当中只考虑位错的等效体力，在计算单元积分时等效体力的计算只有一项：

图 1 将位错源等效为体力示意图图中r点为位错面上任一点，r′为计算区域内任一点，Σ表 示位错所在的几何面，ν为包围位错面的法方向，[ u ]为位错，f 为等效体力.Fig. 1 Sketch of equivalent bodyforcer is a point on the dislocation plane and r′ is a point in the computing domain. Σ is the dislocation plane. ν is normal direction of the interior surface surrounding Σ. [ u ] is the dislocation. f is equivalent force.

Okada的理论(1985，1992)被人们广泛应用于计算地震位错造成的同震形变.这里我们提供一个算例，计算半无限空间中位错源引起的位移以及应力，并将结果与Okada理论的计算结果做对比，以验证我们提出的这种处理位错的方法是否有效，以及计算程序是否编写无误.

如图 2所示，我们计算一个竖直右旋断层所造成的同震变形以及应力.但是Okada的理论基于半无限空间，在有限元的计算当中无法实现.所以这里将计算区域取的足够大用来逼近半无限空间，同时在四周以及底面，使用Okada理论的位移值作为第一类边界条件.模型的上表面使用自由边界条件.

图 2 算例模型示意图 断层为长50 km宽10 km的竖直断层，沿X轴方向位错量为 3 m.计算区域取为1000 km×1000 km×500 km.为了使模拟结果尽量准确，这里将计算区域取的足够大以逼近半无限空间.Fig. 2 Settings of the test case The computing domain is a 1000 km×1000 km×500 km box. A dextral strike slip fault is located at the center of the upper surface. The dimension of the fault is 50 km in horizontal and 10 km in vertical. Slip is in the X direction with a magnitude of 3 m.

图 3显示该算例中与位错等效的体力分布情况，这里画出的是(7)式的积分结果，即对计算区域离散化之后的节点力.等效体力含有δ函数，(7)式的积分只在包含断层的单元内有非零值，因此图中箭头分布在断层周围，在断层两侧箭头的方向与位错方向一致，在断层两端与位错方向垂直.等效体力的这种分布形状符合点位错的双力偶模型，将断层面看成许多点位错组成，最终形成的体力等效于许多双力偶的叠加，于是形成图 3中的形状.

图 3 等效体力分布图 图中箭头表示等效体力的大小和方向，色标表示等效体力的x分量.绿色平面为断层面.Fig. 3 Distribution of equivalent bodyforce in the test case Arrows indicate magnitude and direction of equivalent bodyforce. Color scale indicates the x component of equivalentbodyforce. Green rectangular is the fault plane.

图 4表明，位移计算结果与解析解之间很好吻合，在本文例题采用约40万单元计算时，位移结果的相对误差在大部分区域低于3%，在位移为0的区域相对误差较大.这是因为解析解的结果是0，但是数值解结果即使是一个接近于0的数值波动，也一般不严格等于0.同样的，应力的计算结果也在0值附近有较大的相对误差，但是大部分区域相对误差不大于10%.由于计算应力时需要对位移结果求导，这样会损失一阶精度，所以应力计算的误差要大于位移结果.

图 4 计算结果与解析解对比 (a)为Okada解析解的结果;(b)为数值计算的结果.Fig. 4 Comparison of our numerical result with Okada′s solution (a) and (b)are Okada′s analytical solution and our numerical result respectively.The x component of displacement is showed here.

另外需要指出，断层处的单元尺寸对计算结果会造成一定的影响.这包括两方面的内容：其一，断层处位移梯度大，根据有限元的一般理论，在解的梯度较大的区域，需要使用较精细的网格；其二，虽然等效体力与位错造成的结果完全等价，但是用到有限元方法上来，需要保证离散化之后的节点力尽可能地反映出等效体力，(7)式表明节点力是数值积分的结果，很显然，积分区域剖分的越细，积分结果越精确，于是离散化之后的节点力越能反映等效体力，结果也就越精确.总之在断层附近的单元越精细，计算结果越准确.本文的两个算例都采用在断层附近自适应加密的方法进行计算.

图 5 华北地区网格覆盖情况 图中显示华北地区网格分辨率信息.由于断层处位移梯度 大，所以越靠近断层加密次数越多，网格分辨率越高.华北地区分辨率普遍在20~40 km.Fig. 5 The mesh scales in North China The mesh is densified where the solution has steep gradients. As a result，our mesh is denser near the fault and coarser far from it. The scale of our mesh is between 20 km and 40 km in North China.

在通过初步的测试之后，我们使用该方法计算2011年3月11日发生于日本的Tohoku-Oki特大地震(M9.0)的同震位移及应力，并探讨其对我国华 北地区主要断层上库仑应力的影响.特大地震造成 的位移必须考虑地球的球形效应，因此计算区域为全球，采用PREM模型，边界条件为外表面自由，半径1000 km处固定.使用六面体网格剖分，在计算过程中自适应加密，即在位移梯度大的地方加密网格，以提高计算精度.最终网格规模约440万单元，断层附近分辨率达到约2 km，华北地区网格分辨率约20~40 km.发震断层采用了Simons 等(2011)的位错模型.

图 6显示该算例中的等效体力分布情况.图中画出投影到断层面上的分量，等效体力与断层位错方向一致.但是作为三维矢量，等效体力与断层位错在空间中的方向并不一致，这与前一个算例不同.造成这种现象的原因是断层模型较前一个算例复杂，例如几何形状的不规则以及位错方向并不严格贴合 断层面等.如果使用点位错的双力偶模型来解释，则此例对应三种位错模式的叠加，而体积模量和剪切模量数值上的差异，会使得不同位错模式等效成体力时有不同的比例因子，叠加以后方向就会改变.在弹性力学中，这种施加力的方向与该力造成的位移方向不一致的情形并不罕见，例如单轴拉伸实验，轴向的施力可以导致侧向的收缩.

图 6 等效体力与断层面位错分布图 图中黑色箭头表示离散化之后的等效体力沿断层面的 分量，白色箭头是断层面上盘相对于下盘的位错沿断层面的分量，红-蓝色标表示断层面上位错量的分布.Fig. 6 Distribution of equivalent bodyforce and slip vectors Black and white arrows indicate the projection of equivalent bodyforce and slip vectors onto the fault plane respectively. The color scale indicates the magnitude of slip vectors.

图 7是华北地区位移计算结果与GPS观测对比，方向和大小都有比较好的吻合，进一步证明计算结果的可靠性.

图 7 计算位移与GPS观测位移对比 图中显示华北地区位移计算的结果.黑色箭头表示GPS位移，彩色箭头表示计算位移.Fig. 7 Comparison of our numerical result and GPS observations Black arrows and colorful arrows indicate GPS observations and displacement of our numerical result respectively.

图 8 应力计算结果 这里展示15 km深度处的应力结果，拉张为正.(a)是东西向正应力云图及等值线，(b)是南北向正应力云图及等值线，(c)是剪应力的theta-phi分量(theta与phi分别对应地理坐标系的南、东方向)，(d)图中“十”字型的线段表示水平主应力的方向和大小，曲线是水平最大主应力的等值线.由图可见最大水平主应力为拉张，方向指向断层，大小在千帕量级；最小水平主应力为压缩，量级在百帕以下.Fig. 8 Numerical result of stress fields The stresses at a depth of 15 km are plotted. We adopt the convention in elastic mechanics by which tensile stress is defined to be positive.(a)(b)(c)show the east-west，north-south and theta-phi(theta is corresponding to south and phi is corresponding to east in geography coordinate)component of the stress tensor respectively. In(d)，crosses are plotted according to principal stresses in horizontal plane and curves are contour lines of the maximum component.

由应力场的结果可以进一步计算华北地区主要断层上库仑应力的变化.据库仑应力计算公式ΔCFS=Δ| τ|+μΔ| σn |，其中Δ| τ|、Δ| σn |分别为作用于断层上的剪应力大小和正应力大小的变化量，|σ_n|在拉张情况下取正值，压缩情况下取负值，μ 为内摩擦系数.计算得到的ΔCFS值越大，表示在该 应力场的作用下，断层越趋于失稳，反之则断层越稳 定.库仑应力考虑的只是应力场的增量 Δ|τ|和Δ| σ_n |，不需要知道全部初始应力场.假设华北地区断层的法向在日本地震前后未发生改变，则可直接由Δσ_n计算Δ|σ_n|.但同样的方法不适用于剪应力.为了求出Δ|τ|，我们使用Δτ在初始应力场剪应力方向的投影作为近似，即：Δ|τ|≈Δτcosθ，其中θ为Δτ与初始剪应力方向的夹角，因此需要知道初始应力场的剪应力方向(石耀霖和曹建玲，2010).本文假设华北地区的构造应力场以NEE方向为最大压应力方向，在此基础上计算出各断层面上的剪应力方向，再根据上述讨论计算ΔCFS.计算结果显示3 · 11地震产生的同震应力场，使得我国华北地区大多数 断层趋于稳定，库仑应力增量从几十帕到几百帕不等.

图 9 华北地区主要断层上库仑应力变化 (a)为各断层分布图以及库仑应力云图；(b)为各断层库仑应力数值，其中每条断层上给出不同位置库仑应力变化的最大值和最小值.各断层分别为： 1乌拉山北，2榆林—府谷，3山西断陷，4五台山前，5蔚县—延庆，6尚义—怀安，7张家口—北票西，8张家口—北票中，9张家口—北票东，10紫荆关，11太行山前北，12沧东，13宝坻，14唐山—宁河，15益都，16郯庐.Fig. 9 Coulomb failure stress changes on various faults in North China (a)shows positions of various faults in North China and contours of Coulomb stress changes on them.(b)shows the maximum and minimum of Coulombs stress changes on the faults. The faults are: 1 Northern Ural Mountains，2 Yulin-Fugu，3 Shanxi fault depression zone，4 Front of Wutai mountain，5 Weixian-Yanqing，6 Shangyi-Huai′an，7 Zhangjiakou-Beipiao west，8 Zhangjiakou- Beipiao middle，9 Zhangjiakou-Beipiao east，10 Zijingguan，11 Front and north of Taihang mountain，12 Cangdong，13 Baodi，14 Tangshan-Ninghe，15 Yidu，16 Tanlu.

本文介绍了一种模拟弹性位错问题的等效体力方法，将解不连续的问题转化为解连续的问题.使得建模时可以完全不用理会不连续面，大大降低了建模的难度；由于不用考虑贴合不连续面，也使得单元质量得到一定程度的保证，从而保证求解的精度；另外，由于不同问题的不连续面位置、形状等各有不同，因此含有不连续面的模型不具备可重用性，而使用等效体力方法可以避免含有不连续面，因此建立的模型具有很好的可重用性，一个模型可以处理多个问题.该方法对位错的处理体现在单元积分中，将位错等效为体力项，使用本文公式(4)计算出这个体力的值，然后将此体力项的贡献添加到方程的右端项中去，常用的有限元软件均支持体力项的添加，因此本文介绍的方法很容易实现.

使用该方法计算半无限空间竖直断层的同震位移和同震应力场，并将结果与Okada的解析解结果 对比，位移结果误差约为3%，应力结果误差约为10%.

最后计算了2011日本东北大地震在我国华北 地区产生的同震位移和同震应力.同震位移的计算 结果同GPS观测值符合较好.应力计算结果表明此次地震在我国华北地区造成的主张应力近东西方向，大小达到千帕量级；主压应力接近南北方向，大小达到百帕量级.假设华北地区的构造应力场为NEE方向的挤压，该地区的断层如果发生错动，均为此构造应力场控制的错动.在这个假设之下，我们计算了华北地区各大断层上库仑应力的增量，结果显示各大断层上的库仑应力基本上呈减小态势，量级从几十帕到几百帕不等，这意味着2011日本东北大地震的同震应力场，使我国华北地区各大断层趋于稳定.