地球物理学报  2013, Vol. 56 Issue (6): 2133-2139   PDF    
关于“用单元降刚法探索中国大陆强震远距离跳迁及主体活动区域转移”的讨论——横向各向同性“杀伤单元”才是更好的途径
董培育1,2 , 石耀霖1,2     
1. 中国科学院计算地球动力学重点实验室,北京 100049;
2. 中国科学院大学地球科学学院,北京 100049
摘要: 《地球物理学报》第55卷第1期的"用单元降刚法探索中国大陆强震远距离跳迁及主体活动区域转移"一文提出用减小剪切模量的方法模拟地震断层错动的效应, 这种处理虽然可以模拟压扭性走滑断层错动时的剪应力降低, 但会导致垂直于断层的正应力也剧烈变化, 因此这种简单的减小弹性模量的方法并不合理.本文探讨了一种更好的方法——横向各向同性"杀伤单元", 利用该方法模拟断层滑动效应得到了较为合理的结果, 与Okada解析解吻合较好, 能够正确反映断层滑动时的应力变化.
关键词: 有限元      降刚      剪切模量      横向各向同性      杀伤单元     
A discussion on "The mechanism of long-distance jumping and the migration of main active areas for strong earthquakes occurred in the Chinese continent"——transverse isotropic "wounded element" is a better method
DONG Pei-Yu1,2, SHI Yao-Lin1,2     
1. Key Laboratory of Computational Geodynamics, Chinese Academy of Sciences, Beijing 100049, China;
2. College of Earth Science.University of Chinese Academy of Sciences, Beiing 100049, China
Abstract: In a paper titled "The mechanism of long-distance jumping and the migration of main active areas for strong earthquakes occurred in the Chinese continent "(Vol.55, No.1 of Chinese J.Geophys.(Chinese Edition)), Yang et a1.proposed to simulate the effects of earthquake fault slip by reducing shear modulus of the elements in calculation.Although this method can simulate the shear stress drop of a fault under compression and shear, it will produce a significant drop of normal stress on the fault which is incorrect.In this paper, we propose a better method——transverse isotropic "wounded element", and get a reasonable result, which is in good agreement with analytic solutions in the stress calculation as a result of fault slip..
Key words: Finite element method      Reduced element stiffness      Shear modulus      Transverse isotropic      Wounded element     

在《地球物理学报》第55卷第1期的“用单元降刚法探索中国大陆强震远距离跳迁及主体活动区域转移”[1]一文中提出:震源区部分丧失承载能力是引起强震远距离跳迁及主体活动区域转移的重要原因,而在模拟计算中采用显著降低模型中当前包含最大应力单元的一组单元的弹性模量的方法,计算模拟强震断层活动引起的应力场的变化调整.调整后的应力场中再次寻找包含最高应力单元的单元组,作为后续强震的发震位置,再次显著降低这组单元的弹性模量,又引起新的应力应变场的调整和应力集中转移.与此同时,上一次被降低了弹性模量的单元组恢复性增加其弹性模量数值,模拟断层的愈合.

在有限元计算中,研究张裂隙破裂,或工程中的开挖造成的应力变化,利用“杀死单元”(kill elemem)——即假定裂开或挖走的单元杨氏模量降为0,是一种得到广泛应用的方法[2-3]。但地震剪切破裂时,虽然有应力降,但应力并不降为0,断层仍具有一定的抗剪强度,因此如何将杀死单元的思路应用于地震断层,是一个有意义的课题,文献[1]的这种方法是一次非常有益的尝试探索.通过文献[1]“杀伤”而非“杀死”的方法,实现使震源区域部分丧失承载能力,产生应力降[1, 4-5].但是该文这种简单的减小弹性模量的方法也还存在严重的问题.根据各向同性下的虎克定律,如果泊松比ν不变,弹性模量E和剪切模量G是相互联系的,降低了断层带内的弹性模量E,剪切模量G=E/(2(1+ν))也会同时降低;或者说降低了剪切模量G,弹性模量E也会按比例降低.用降低剪切模量G模拟断层剪切破裂是一种看似可行的思路,但同时弹性模量E降低,这不能正确模拟断层在正应力下的作用.文献[1]中提到“单元的泊松比也做相应的变化,单元的剪切模量和体积模量合理.E的原值是80 GPa,减弱到17.8 GPa”.如果令G保持不变,ν=E/(2G)-1,那么大幅度减弱E之后的ν变为负值,虽然聚合物泡沫材料发现有泊松比可以是负值的材料[6],但是在自然条件下的岩石介质中负值的泊松比基本没有.所以即便是泊松比取值在合理范围内相应调整,总是会在减弱E的同时,G也减弱,或减弱G同时,E也减弱.

我们利用一个简单的二维的平面应力问题的例子来探讨这个问题:假设初始状态没有断层的存在,整个区域都是均匀各向同性的.在边界荷载的作用下,发生剪切破裂,产生走滑断层,如图 1所示,假设断层长1 km,厚度为2 m (灰色区域为断层,断层夸大).边界条件为:左边界。垂直向下的位移2 m,水平向右的位移1 m;右边界。垂直向上的位移2 m,水平向左的位移1 m;上、下边界的条件相同,均为垂直方向自由,水平方向从左到右位移线性变化.这样的边界条件既有左旋剪切的作用,也有东西向垂直于断层的正压应力的作用.

图 1 计算采用的含有断层的模型 Fig. 1 Calculation model which contains fault

把上述边界条件下的均匀各向同性不含断层的二维模型的位移和应力作为断裂前的初始状态,与后面断层发生滑动后的状况比较.考虑断层发生时,断层带内的物质采用文献[1]的方法,减弱其弹性模量E,计算其在上述边界条件下的位移和应力状态(图 2).计算值减去初始状态的值得到的差值,反映的即是发生地震前后的状态变化(图 3),这一差值可与Okata的解析解(图 4)进行比较.

图 2 文献[1]单元降刚以后计算的位移和应力 (a)x方向(垂直于断层)的位移Ux;(b)y方向(平行于断层)位移Uy;(c)剪应力σxy(d)正应力σxx. Fig. 2 The displacement and stress distribution after reducing the Yong'S modulus of fault elements with the method used in Ref. [1] (a) Displacement Ux of the direction x(vertical to the fault); (b) Displacement Uy of the direction y(parallel to the fault); (c) Shear stress σxy; (d) Normal stressσxy.
图 3 文献[1]单元降刚以后计算断层发生后与发生前位移和应力变化量 (a) △Ux;(b)△Uy;(c)△σxy;(d)σxx. Fig. 3 The change of displacement and stress before and after fracture with the method used in Ref.[1]
图 4 Okata解析解计算的地震发生后位移和应力变化量 (a)△Ux;(b)△Uy;(c)△σxy;(d)△σxx. Fig. 4 The change of displacement and stress before and after fracture with the Okada's analytical solution

图 2中虽然反映出断层走滑对y方向位移和剪应力的影响,但是由于断层的杨氏模量被人为地降低,导致x方向位移出现了不合理的状况,断层西面和东面的岩体均向断层移动(图 2a),造成了断层东西部广大区域的东西向压应力下降(图 2d).而走滑断层发生应该主要降低断层的剪应力,不应该对断层上正应力有如此巨大影响的.

这一问题在图 3图 4的比较中看得更为明显.图 3是用文献[1]的降刚法计算的地震前后变化,图 4是在断层在相应位错量下用Okata公式计算的地震后位移和应力变化量.图 4d显示Okada解地震前后正应力σxx几乎没有变化,而图 3dσxx却有兆帕量级的降低.这主要在于文献[1]降刚法使断层东西两侧广大区域在x方向发生了实际不应该有的位移(图 3a),这也导致剪应力降与Okada解有显著差别(图 3c图 4c比较).

我们建议采用横向各向同性的“杀伤单元”(wounded element)来处理断层破裂问题,即把达到破裂的断层单元的力学性质改为弱化了的横向各向同性物质,平行于断层剪切面的剪切模量降低,从而模拟断层发生时的错动和应力降;但垂直于断层面的杨氏模量不变,模拟法向应力基本不变的事实.横向各向同性物质,如图 5所示,假设1轴的方向为对称轴,横向(即2-3平面内)为各向同性的,但是沿着1轴的法向变形与2或3轴方向可以不同(Ei代表沿着i轴的弹性模量),1-2、1-3平面内剪切变形与2-3平面内剪切变形也可以不同(G21G31相同,但与G23不同),νij代表i方向上受力时,与它垂直的j方向上相应变形的泊松比;Gij代表法向量为i的平面上沿着j轴方向的剪切模量.这时共有5个独立参量,本构关系见图 5[7].这里由于我们采用二维平面应力问题以1-2平面作为剖面,则所需的独立弹性参数进一步减小为4个,即E1E2ν21G21.

图 5 横向各向同性材料图示和本构关系方程(注意其中的,因此只有5个独立参量) Fig. 5 The sketch of the transverse isotropic materialand theconstitutive relation in which SO that there is only 5 independent parameter.

在模拟断层发生时,背景单元仍然为各向同性,我们仅仅降低断层单元的剪切模量G21,垂直于断层方向的杨氏模量E1并不需要变化,即E1=E2,均取原来各向同性岩石的E的值,泊松比ν21也取原来各向同性岩石的值ν0.然后进行同样的计算,得到断层的滑动量,如图 6所示,断层发生后应力状态的变化见图 7.把图 7与给定类似断层错动,根据Okada给出的解析解[8],利用汪荣江的开源计算程序[8-9]来计算的图 4比较,可以看出位移和应力变化量,不但空间分布图像相似,而且各点计算值大小也很相近,各相应图件仅有微小差别.应该指出,Okata的解适用于厚度趋于0的线状断层,我们的断层有有限厚度,二者不会完全相同.另外目前我们计算中把区域简单按矩形网格划分,因此对断层尖端附近的应力集中区域,网格不够密,数值计算误差会比较大,如果改进单元划分,采用非结构化网格对断层端点附近加密,应该会取得更佳结果.但目前的计算已经足以说明,我们建议采用的横向各向同性“杀伤单元”减小G21的方法,能够较好地反映断层滑动时的应力变化.

图 6 用横向各向同性介质消减G21模拟断层滑动的计算结果 (a)Ux;(b)Uy;(c)σxy;(d)σxx. Fig. 6 The result of simulation fault slipping with the method of reducing G21 in the transverse isotropic material
图 7 横向各向同性介质消减G21模拟断层滑动的计算结果 (a)△Ux;(b)△Uy;(c)σxy;(d)σxx. Fig. 7 The result of simulation fault slipping with the method of reducing G21 in the transverse isotropic material

总之,强震跳迁过程的地球动力学模拟是非常复杂困难的课题,文献[1]的基本思路是有价值的,采用单元降刚,即显著降低其弹性模量的方法,是一种新的探索,有助于了解大范围的兆帕量级的应力场的调整,解释后续强震远距离跳迁的问题.但是在存在垂直断层的正应力情况下,简单通过减弱其弹性模量E的具体方法处理剪切断层存在严重的问题.我们建议采用的横向各向同性“杀伤单元”减小G21,的方法,能够较好地反映断层滑动时的应力变化.而随时间演变逐步增加G21,则能够反映断层的愈合.

致谢

本文得到程惠红和尹凤玲的有益建议,特此感谢.

参考文献
[1] 杨树新, 陆远忠, 陈连旺, 等. 用单元降刚法探索中国大陆强震远距离跳迁及主体活动区域转移. 地球物理学报 , 2012, 55(1): 1105–116. Yang S X, Lu Y Z, Chen L W, et al. The mechanism of long-distance jumping and the migration of main active areas for strong earthquakes occurred in the Chinese continent. Chinese J.Geophys. (in Chinese) , 2012, 55(1): 1105-116.
[2] 阎波, 方云. 基于ANSYS的边坡开挖模拟. 山西建筑 , 2008, 34(1): 93–94. Yan B, Fang Y. Simulation ot slope excavation based upon ANSYS. Shanxi Architecture (in Chinese) , 2008, 34(1): 93-94.
[3] 谢海国, 章广成, 杨昌斌. 边坡分级开挖过程的动态模拟. 土工基础 , 2007, 21(3): 54–56. Xie H G, Zhang G C, Yang C B. Achievement of dynamic simulation of excavation of slope in a hydropower station. Soil Eng.and Foundation (in Chinese) , 2007, 21(3): 54-56.
[4] Lu Y, Yang S, Chen L, et al. Mechanism of the spatial distribution and migration of the strong earthquakes in China inferred from numerical simulation. Journal of Asian Earth Sciences , 2011, 40(4): 990-1001. DOI:10.1016/j.jseaes.2010.07.008
[5] 陆远忠, 叶金铎, 蒋淳, 等. 中国强震前兆地震活动图像机理的三维数值模拟研究. 地球物理学报 , 2007, 50(2): 499–508. Lu Y Z, Ye J D, Jiang C. 3D numefial simulation on the mechanism or precursory seismicity pattern before strong earthquake in China. Chinese j.Geophys (in Chinese) , 2007, 50(2): 499-508.
[6] Lakes R S. Foam structures with S negative Poisson'S ratio. Science , 1987, 235(4792): 1038-1040. DOI:10.1126/science.235.4792.1038
[7] Slawinski M A. Waves and Rays in Elastic Continua. (2nd Ed). World Scientific, 2007 .
[8] Okada Y. Surface deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America , 1985, 75(4): 1135-1154.
[9] Wang R, Martin F L., Roth F. Computation of deformation induced by earthquakes in a multi-, layered elastic crust——FORTRAN programs EDGRN/EDCMP. Computers & Geosciences , 2002, 29(2): 195-207.