有限元法主要原理为用有限个单元所构成的离散化结构代替原来的连续体结构,来分析构造的应力与变形。与其它方法相比,有限元法具有适用性强的特点,即无论待计算的材料性质与外载荷如何变化,均可应用该方法进行计算(李人宪,2002; 龙志飞,岑松,2001)。

(1)一般性研究方法

有限元法应用于构造应力场与地震预测研究中,其一般性的研究方法与步骤为:简化所研究目标区域的活动断层; 建立一个符合实际的二维或三维有限元模型; 施加合理的模型边界条件(位移与应力等); 利用有限元计算求解出目标区的应力应变场、应变能密度等参数; 用有限元后处理获取位移、应力场结果,并可进一步实现结果可视化; 对数值模拟计算结果进行综合研究分析。

(2)主要处理环节

将有限元法应用于构造应力场与地震预测研究的过程中,需要考虑的主要处理环节有:断层处理方法、介质物性参数的选取、边界条件的加载和数值积分求解方法等。

① 断层处理方法

基于有限元理论方法与计算条件的限制,以往大多数方法是采用连续体模型进行分析的,仅将断层处理为软弱地带。但是随着不连续体理论与计算算法的发展(余天堂,2007),不连续模型逐渐开始应用于构造应力场与地震预测研究中。依据断层处理方法不同,模型建立可分为两类:(a)将研究区作为一个连续体,仅将断层处理成软弱地带,并将不同地块与断层软弱地带赋予不同弹性模量与泊松比,不考虑断层间的作用;(b)将研究区作为一个不连续体,考虑断层因素,采用非连续接触分析来处理断层间的非线性接触关系。

② 介质物性参数选取

根据所研究问题的不同,将研究区地质模型进行横向分区与纵向分块,根据实际情况,将不同地块考虑为弹性、粘性、粘弹性、弹塑性等不同介质,并选取不同的物性参数(傅容珊,黄建华,2001),具体参数如表1所示。

表1 常用介质物性参数对比
Tab.1 Comparison of petrophysical parameters of common media

③ 边界条件加载

有限元计算的边界条件主要采取以下两种方式进行加载:(a)利用研究区的GPS观测数据,通过插值计算获取边界GPS数据,把在研究区域的边界加载的边界位移或位移速度,作为边界条件;(b)根据参考震源机制解与地表应力测量结果,采用试错法反复调整模型边界作用载荷的大小和方向,再与获得的研究区域GPS结果进行比较,尽量使模拟得到的结果与最新获得的GPS观测资料基本相吻合,最终确定出模型边界的加载条件。

④ 数值积分求解方法

数值积分法是指不通过坐标变换,直接用数值方法对动力平衡方程积分。其求解方法主要分为显式算法和隐式算法两类(刘立终,刘相华,2001; 李初晔,王增新,2010)。显式算法基于动力学方程,不需要进行平衡迭代,计算速度快,时间步长只要取的足够小即可,一般不存在收敛性问题。故其所需内存比隐式算法少,且具有较好的稳定性; 而隐式算法基于虚功原理,需要迭代过程,其每一增量步内都需要对静态平衡方程进行迭代求解,所需计算内存较大,且易受到迭代次数及非线性程度的限制,计算效率不高。

20世纪70年代中期,我国学者将有限元法引入构造应力分析中(罗焕炎,1974),并详细阐述了有限元法在地质力学分析中的应用。随后,有限元法在地球科学定量分析研究中迅速发展,应用越来越广泛。现将有限元法在构造应力场与变形场、地震应力场与地震预测等方面的研究历史与现状分述如下。

(1)构造应力场与变形场的模拟

汪素云等(1980)对中国及邻区现代构造应力场进行了数值模拟,通过五种应力和位移边界条件的试算,获得合理参数,得到了我国构造应力场分布状态。许寿椿和朱正(1983)采用二维非线性有限元法,根据郯庐断裂区的原地应力测量资料反演了区域应力场方向。杨光宇(1985)采用平面问题的有限元法,对川滇地区强震活动与构造应力场的关系作了初步研究。梁海华(1987)运用有限元法模拟汾渭断陷带构造的地质和地球物理特征,并探讨了其动力条件和构造应力分布。宋惠珍等(1987)对北京地区区域应力场进行了数值模拟,其结果为区域应力场提供新的资料和证据。强祖基和谢富仁(1988)根据断层错动资料和小地震的平均震源机制解结果,运用弹性力学的叠加原理,对临汾裂谷作三维有限元模拟,反演其构造应力场,结果表明,临汾裂谷现今主要受南东东—北西西向引张力作用,辅之以北东—南西方向的挤压,同时在裂谷的中心叠加有底层拖曳的深部作用。张东宁和高龙生(1989)利用三维弹性有限元方法对欧亚板块东部进行了数值模拟,计算结果认为:水平向的最大主压应力轴走向在中国大陆地区呈辐射状分布,其中心大体上在青藏高原。张东宁和许忠淮(1994)对青藏高原进行三维弹—粘性有限单元模拟,获得高原岩石圈的构造应力状态及其物质的运动特征。模拟结果显示青藏高原的重力势能是中国大陆现代构造特征的主动力。丁原章和黄新辉(1995)采用二维有限元分析方法反演南海北部及其邻近地区的构造应力场,探讨可能的发震危险地段,为地震工程提供地质背景资料。徐志斌等(1998)对晋中南地区的构造应力场进行数值模拟,结果表明,现代构造应力场差应力在沁水含煤向斜轴部存在着呈北东向展布的高值条带,且最大主应力为压应力。汪素云等(1996)利用中国及其邻区地震震源机制解推断的观测应力方向,反演了周围板块作用力的相对大小。武红岭等(1996)采用粘性牛顿流体有限元数值方法对郯庐断裂带主要构造演化阶段的应力场,形变场和运动速率进行数值模拟,讨论了各构造时期的动力来源、边界力的作用方式及岩石物性影响。尹京苑和梅世蓉(1998)总结了三维数值模拟计算中在边界条件选择上存在的若干问题,得到板内地壳中3个主应力和3个最大剪应力随深度的变化曲线。

刘科等(2002)采用平面应力的有限元法反演河北平原区的构造应力场,确定最大主应力、最大剪应力集中区。宋惠珍等(1999)采用二维显式有限差分法与三维隐式有限元法,把断层面处理成非连续面,分别计算了模型剖面上断层逆冲运动和三维断层逆冲运动伴随的断层邻域应力场。Cai等(2000)提出了LDDA(Lagrangian Dis-continous Deformation Analysis)方法,对块体之间接触采用了摩擦接触准则,这种方法可动态模拟出地震过程。傅容珊等(2002)研究了印度板块与欧亚板块碰撞40 Ma以来青藏高原挤压隆升、东亚大陆形变以及构造应力场的演化过程。模拟所得水平形变速度与GPS观测的格局很吻合。肖兰喜等(2003)给出印度板块、太平洋板块和菲律宾板块对中国大陆驱动的边界作用强度之比约是4:1.25:1。曹雪峰等(2003)运用有限元法,对粘弹介质状态下的南北地震带及其相邻区域应力场时空演变过程作了数值模拟研究。计算出了十万零十年累积载荷作用下该模型应力场的演变。李延兴等(2004)利用推导出的块体弹性运动方程与块体GPS站速度,计算并分析了中国大陆及周边地区应变场的基本特征,所得主压应变方向与用地质学方法和测震学方法得到的主压应力轴方向具有很好的一致性。郑文衡等(2005)则采用非线性动力学模型模拟地震活动,初步研究了地震动态触发机制。陈连旺等(2005)建立了中国大陆构造应力应变场的二维有限元模型,通过数值模拟方法研究了中国大陆现今构造应力应变场年变化图像的基本特征。陈晓利等(2005)采用二维有限元法模拟了渤海海域现代构造应力场,得出其各个区域的应力场特征。王连捷等(2006)对青藏高原中段进行了构造应力应变场的数值模拟,结果显示,其速度场由南向北逐渐减小,且两侧逐渐减小,而最大主应力总体近SN向。章纯(2007)运用有限元数值模拟方法,结合板块构造边界条件和板内构造分布特征对中国大陆东部地区的基本构造应力场进行数值模拟,分析板块边界作用力的变化对内部应力场的影响。李细光等(2007)对广西及其邻区的构造应力场进行了有限元模拟,通过该区域的应力应变等分布特征,得出广西地区构造应力场总体上具有分区性,桂东南地区为高应力集中区,桂中地区为弱应力分布区。

(2)有限元法在地震应力场与地震预测中的模拟

1976年唐山地震后,为了预报未来地震,王仁等(1980)首先进行了相关数值计算,通过对1966年邢台地震以来华北地区的5次大地震后构造应力场变化的模拟,预测未来的地震危险区。这个方法曾被用于首都圈地区(王仁等,1982; 宋惠珍等,1982)、鲁南地区(张郢珍等,1988)、云南地区(郑治真等,1984)、西南地区(杨光宇,1981; 杨光宇,1982; 王继存等,1991; 张周术等,1991)的地震危险区模拟。其中,数值计算均采用弹塑性模型,初步推测出震后应力场危险区的范围和可能震级的大小。

朱岳清等(1988)提出了唐山地震孕震的物理力学模型,采用有限元法对孕震过程进行三维力学分析,并认为摩擦滑动破坏机制比较符合唐山地震的实际情况。郑治真等(1988)采用有限元法模拟在压应力作用下的层理分布和应力集中区,其数学模拟的应力集中区与台湾省的强地震分布区基本一致。陶振宇和张立明(1991)根据构造地震的发生机制建立了一个数学模型。该模型不仅能模拟岩石摩擦试验的粘滑特性,而且在该模型基础上建立起来的数值模拟方法能对构造地震的发生作出定性和定量的分析,得到了一些对震源深部岩石力学性质的认识。

滕春凯等(1992)使用有限元法研究了含摩擦多断层问题,结果表明,断层内的摩擦力对应力场的分布有很大的影响,并可能对断层的扩展行为有着非常重要的影响。丁原章等(1992)采用三维弹性、粘弹性程序,模拟广东新丰江水库区地壳上、下层的构造骨架,探讨断裂的危险程度和下层构造变动对上层构造以及地形地貌的影响。孙荀英等(1994)采用三维粘弹性有限元方法拟合唐山地震震时和震后的水平与垂直地形变,反演华北板块下方深部物质的流变学性质。Shi和Zhang(1994),耿鲁明等(1994)专家进一步应用非线性动力学模型模拟了地震活动的研究。

王妙月等(1999)从线性流变体介质内制约质点运动的运动方程出发,导出了模拟地震孕育、发展、发生动态过程的三维有限元方程及程序,并完整地模拟地震孕育、发展、发生的全过程,为地震过程本质的认识以及物理预报的实现提供了一个潜在的新手段。Cai等(2000)采用LDDA有限元法成功模拟了唐山地震动力学过程,但是只有降低摩擦系数才能产生地震。朱守彪和石耀霖(2004)根据震源机制解和地质调查资料,采用伪三维遗传有限元法反演了中国川滇部分地区受到的边界作用和该地区底部所受的剪切作用力。结果认为,青藏高原物质受挤压向东和东南运动过程中,下地壳物质比上地壳更易于流动,从而对川滇菱形块体上地壳有拖曳作用。

2008年5月12日,在四川省龙门山断裂带发生了汶川M_S8.0强烈破坏性大地震。在历史地震活动相对较弱的龙门山地区突然发生了如此强烈的地震,其动力成因是什么?如何累积如此大的能量?震后,众多学者赴灾区进行了科学考察,取得了丰富的第一手资料和研究成果(徐锡伟等,2008; 张培震等,2008; 冉勇康等,2008; 刘启元等,2009; 陈桂华等,2009)。这些资料及成果为深入研究汶川地震的动力学机制奠定了强有力的基础。目前许多专家已初步构建了多种模型来解释汶川地震的发震构造特征,采用有限元等数值模拟方法深入研究并解释汶川地震成因。

王卫民等(2008)根据地质资料和地震形成的地表破裂轨迹,构造了一个双“铲状”有限地震断层模型,利用反演技术重建地震的破裂过程,重现汶川大地震发生时地震断层的滑动情况。其结果显示,汶川大地震主要是沿龙门山构造带的映秀—北川断裂和灌县—江油断裂发生的逆冲兼右旋走滑破裂事件,层面滑动分布显示两个高滑动区先后发生在地震破坏最为严重的映秀和北川地区,最大滑动量高达1 200～1 250 cm,且破裂过程也显示了一定的复杂性。地震破裂的平均走滑量略大于平均倾滑量,与多种观测资料获得的震前龙门山断裂带构造变形相一致(张培震等,2008),推断是由于长期区域应力场作用和龙门山地区特殊的物质组成和结构孕育了这次千年尺度的强烈地震。

陈连旺等(2008)充分考虑了对地质构造运动和地震活动起重要作用的活动断裂,采用接触摩擦分析理论处理活动断裂的力学特征,建立了川滇地区活动断裂三维非线性有限元模型。利用陈运泰等给出的汶川地震主震破裂过程的研究成果,将主震破裂事件分成7个破裂子事件,对应数值模拟计算的7个载荷步。在此基础上,对于汶川M_S8.0地震主震破裂引起的应力应变场变化的空间分布特征及其控制因素,开展了三维非线性有限元数值模拟分析。研究得到汶川地震余震序列空间分布特征及其震源机制类型的力学解释; 汶川地震发生后的龙门山断裂带应力状态的变化特征; 汶川地震导致川滇地区主要活动断裂构造运动的加载效应等。

朱守彪和张培震(2009)采用粘弹性接触问题的有限元法,并考虑重力作用,对青藏高原东缘的应力场空间分布及其随时间的演化进行了数值模拟,结果显示应力在空间由分散分布逐渐向龙门山及周边地区转移集中。他们基于前人的研究成果及计算分析,初步分析了汶川地震孕育发生的动力学过程:青藏高原的东流物质在向东运动过程中由于受到稳定的四川盆地的阻挡,一部分东流物质在川西地区囤积,造成龙门山隆升; 高角度(50°～70°)、犁状的龙门山断层面上的正应力随着川西高原向东运动而不断增大,导致该断层的闭锁性逐步加强,并且分布在断层附近的变质杂岩为存贮高密度弹性应变能提供物质保障。但另一方面随着青藏高原较柔软的下地壳物质的不断向东运动,囤积的东流物质对龙门山断裂带上盘的推挤作用不断加强,从而导致断裂带上剪应力越来越大; 当剪应力超过摩擦强度时,断层解锁产生滑动,发生地震。其模拟结果还表明龙门山断层面上的摩擦系数较高,断裂带上地震的平均复发周期约为3 163 a,这与其他学者的研究结果(张培震等,2008)有一致性。

陈祖安等(2009)用三维流变非连续变形与有限元相结合(DDA+FEM)的方法,在青藏高原及邻区三维构造块体相互制约的大背景中,考虑了龙门山断裂带东西两侧地势、地壳厚度和分层的明显变化,及断裂带东侧四川盆地及鄂尔多斯块体坚硬地壳阻挡的影响,通过用GPS资料做位移速率边界约束和震源机制约束,计算得到研究区的速度场和应力场与该地区GPS结果和震源机制分布结果基本一致。在此基础上,模拟计算现今构造块体边界断层上表征剪应力及法向应力等综合影响的危险度分布,结果表明,上、中地壳层危险度分布中,危险度较高的地段多数与近几十年来发生的七级以上大震区域基本一致。包括2008年汶川8.0级等大震的发震断层。其数值模拟结果显示了汶川大地震的孕育发生机理,并为青藏川滇地区的地震危险性提供参考依据。

汶川地震发生后,众多学者对汶川地震的成因、孕育发生过程及其发生对周边断层地震活动的影响作了大量野外考察与室内数值模拟分析,但上述研究仍处于初级阶段,只有对研究区做进一步的地质、地球化学详细震学反演和高分辨率地球物理勘探等基础工作,并客观地构造物理建模和数学力学等综合分析,才有可能深入了解汶川地震破裂过程,深化认识发震构造模型、发震断层长期滑动习性和地震成因。