本发明公开了一种火山和地震粘弹性形变自动建模有限元计算方法，利用等效体力公式将地震剪切位错和火山膨胀压力源等效为有限元计算中的体力项，再通过局部网格自适应加密技术在地震与火山源处产生足够密的网格，即可省略传统地震和火山形变有限元计算中最费时的前处理步骤‑‑划分断层/火山岩浆房网格，实现火山和地震粘弹性形变的自动建模功能。本发明的有益效果是：该发明一种火山和地震粘弹性形变自动建模有限元计算方法自动完成有限元前处理和计算，计算方法高度并行化，通过软件实现自动处理，既节省建模时间又保证计算精度。
1.一种火山和地震粘弹性形变自动建模有限元计算方法，其特征在于：利用等效体力公式将地震剪切位错和火山膨胀压力源等效为有限元计算中的体力项，再通过局部网格自适应加密技术在地震与火山源处产生足够密的网格，省略传统地震和火山形变有限元计算中最费时的前处理步骤，自动完成建模前处理到有限元计算的全部步骤，该方法具体步骤如下：步骤一：地震剪切位错与火山膨胀压力源的等效体力替代
Mogi模型代表的球状腔体压力源问题用3个正交扩张点源或3个正交张拉位错，这些概念上不同的源可以产生源外部相同的位移场；3对正交力偶Fh产生的位移为：ui(x)＝-MjkGik,j(x) [1]其中Gik(x)是格林函数，代表源处k方向单位力产生的点接收点x处的i分量位移，并且矩张量为Mjk＝Fhδjk；代入弹性半空间内部集中力的格林函数表达式，当岩浆房半径a和岩浆房压力变化ΔP满足Mogi模型假设，则存在3个正交力偶，它们产生的位移场与Mogi模型的位移场等效；且由3个正交力偶代表的体力项，可以写成：其中，δx＝x'是狄拉克函数，代表x'处(源处)的点源集中力，▽δx＝x'代表一对符号相反的脉冲，每一对力偶的强度可表示为：由于狄拉克函数可用高斯分布函数的极限来近似，因此将体力项表达为高斯分布函数的形式如下：其中，αxi是高斯函数在xi方向的变量，αxi取值主要取决于网格大小；式[4]中的αx1αx2αx3π3/2代表了单位体积内高斯函数的正则化系数，这种处理保证了等效体力函数足够光滑；
并且，体力函数可以用来模拟不同的岩浆房形状，当αx1＝a,αx2＝b,αx3＝c且α≠b≠c时，式[5]可以用来模拟椭球状火山膨胀压力源，且椭球半轴分别位a，b和c；其他形状的压力源，可以根据弹性力学叠加原理，用式[5]的权重组合形式，即地震矩张量来表示：地震剪切位错的等效体力替代如下：
对于与地震错动等效的剪切位错，与其等效的体力项如下：
其中， 为等效体力，Φm表示单元K中第m个形函数，vj是位错面法向，u(r)i为位错，Cijpq为弹性系数，dΣ是单元内部位错面；
网格自适应加密后，会引入悬挂节点；对网格自适应加密之后悬点的处理方法如下：为了满足解的连续性，悬挂节点8和9的解需要满足式[7]：
上述u8和u9需要从有限元线性系统里消去，待求解完线性系统后通过回带获得它们的值；悬挂点的约束方程可以写成式[8]：uiconstrained＝Cijuj [8]其中，
此时，有限元线性系统，Au＝b [10]需要转换为：constrained
ui ＝Cijuj [12]待求解式[11]后，代入式[12]得到悬挂节点的解；
在计算模型中实现地形起伏和介质横向不均匀性采用不同的方法实现；对于地形起伏，首先生成无地形的球形地球，然后根据高程数据，计算出地表每个节点移动位移量，将目标节点移动到目标位置；为避免移动节点造成网格畸变，同时对地球内部节点施加移动量，内部节点移动量的确定是具体方法是将地表节点位移量作为为地球模型Poisson方程的边界条件，见式[13]，得到地球内部各节点的移动量，Poisson方程解的光滑性保证了使用这种方式添加地形不会造成太大的网格畸变；
对横向不均匀性以及Moho面起伏的处理，不采用移动网格的方式，而是在组建刚度矩阵时，根据Gauss积分点坐标插值得到各积分点的拉美系数和粘滞系数；
利用玻尔兹曼叠加原理，得到maxwell体的积分型本构关系，再利用小应变假设，推导得到小应变假设下的应力递推公式：σ(tn)＝E'ε'+σ0'
σ0'＝β(Δt)σd(tn-1)+σV(tn-1)。
2.根据权利要求1所述的火山和地震粘弹性形变自动建模有限元计算方法，其特征在于：对火山膨胀压力源的体力等效，通过正交力偶和狄拉克函数的引入，得到体力项的力学表达，进一步利用高斯函数对狄拉克函数的等效，克服狄拉克函数难以在有限元方法中处理的难题，用高斯分布函数和地震矩张量的乘积得到火山膨胀压力源的体力等效项；高斯分布函数只与火山膨胀压力源的位置坐标和岩浆房三个主轴有关；地震矩张量的确定需要利用对角单位矩阵进行主轴缩放得到主轴和笛卡尔坐标系一致的椭球状火山膨胀压力源的张量表示，再进行坐标旋转得到任意主轴方向椭球状火山膨胀压力源的地震矩张量。
3.根据权利要求1所述的火山和地震粘弹性形变自动建模有限元计算方法，其特征在于：对地震剪切位错的体力等效，先计算六面体单元被平面切割后得到的多边形，在此基础上计算地震剪切位错在多边形上的面积分，得到地震剪切位错源的等效体力项。
4.根据权利要求1所述的火山和地震粘弹性形变自动建模有限元计算方法，其特征在于：网格自适应加密技术，采用悬挂节点的方式，这种处理方式是一种局部加密，不影响网格的整体复杂度；悬挂节点作为单元新引入的中节点，为了保证解的连续性，在有限元线性系统里面增加额外的约束方程来约束悬挂点；自适应加密进行多次迭代，但每个单元的每个面最多被细分一次，以此来保证不增加网格复杂度，即母单元细分成子单元后，新一轮的加密只能对子单元施加。
5.根据权利要求1所述的火山和地震粘弹性形变自动建模有限元计算方法，其特征在于：基于消息传递界面的高度并行计算方法，先用主进程生成网格，然后采用METIS开源库对网格进行并行分区，对分区后的网格在各子进程上进行刚度矩阵的并行组装，然后通过进程间的通讯组装整体刚度矩阵，使用开源线性方程求解库Trillinos进行大型稀疏线性方程组的并行求解，最后各子进程将各自的计算结果进行结果文件写操作，实现后处理结果的并行显示。
火山和地震粘弹性形变自动建模有限元计算方法
技术领域
[0001] 本发明涉及地球物理领域，具体地说，涉及一种火山和地震粘弹性形变自动建模有限元计算方法。主要应用于测震学和火山形变学领域，其功能是在无需划出断层和岩浆房网格的情况下，利用有限元计算因为地震静态滑动分布/火山膨胀压力源引起的形变以及应力变化，完全实现自动化，节省大量前处理建模时间。
背景技术
[0002] 地震、火山是地球科学的重要研究对象，二者带来的巨大破坏力都对人类生存具有潜在威胁，但目前相关的解析计算方法仍存在较大局限性：大地震往往发生于俯冲带地区或者板内的地形和介质突变区，比如2010年Maule地震、1960年智利Valdivia MW9.5地震、2004年Sumatra Mw9.3地震、2011年的Tohoku-Oki Mw9.0地震和2008年汶川地震等。这些区域往往地形剧烈起伏、陆壳和洋壳厚度差异巨大、介质性质具有很大的非均匀性，如果使用传统的解析/半解析位错理论，如：Okada弹性半无限空间解析式(Okada,1985；1992)和球形位错理论(Sun and Okubo,1998)将无法考虑地表地形起伏、Moho面等地球重要圈层的起伏以及介质的强非均匀性。此外，目前对于火山区岩浆囊压力源引起的地表变形研究多采用解析的Mogi模型，该方法具有表达式简单、计算快速的优点，但Mogi模型无法考虑地表地形，而真实火山区往往具有较大地表地形起伏，例如长白山火山区最高峰海拔约2691m，日本富士山海拔为3776米。综上所述，采用数值方法能够克服解析方法带来的局限性，考虑地表地形、Moho面起伏、介质非均匀性，在地震和火山应变、应力计算中尽可能模拟真实地球，对地震危险性分析和火山区岩浆活动的认知具有重要意义。
[0003] 但目前通用的地震和火山形变的有限元计算方法仍然存在改进空间，最突出的问题就是前处理网格划分十分繁琐。在传统有限元方法里，地震破裂的模拟往往需要在前处理中划出具体的断层破裂面，而火山地区的地下压力膨胀源一般要在前处理中刻画出球/椭球的岩浆房网格然后施加面力边界条件。但这种前处理方式存在以下困难：1、当断层数量较少、且断层主要为高倾角或者近直立时，前处理建模难度相对较低，但当断层倾角较低或者断层数量较多时，前处理往往十分费时且网格质量不能保证；2、当岩浆房形状复杂、或者岩浆房数量较多时，前处理建模难度大，且往往只能采用四面体网格，计算精度低于六面体网格。
[0004] 鉴于以上实际应用原因，非常有必要发展一种新的火山和地震粘弹性形变有限元计算方法，要求这种计算方法自动完成有限元前处理和计算，既节省建模时间又保证计算精度。
[0005] 如果还按照传统有限元方法，寻找一种策略能让程序自动划分出断层和岩浆房网格，显然这种策略的难度很高，而且难以保证网格不过度畸形。因此，我们采用另一种策略，即把断层(剪切位错)和岩浆房膨胀(膨胀位错)等效为体力。此方法不需要像传统有限元设置内边界条件(对于地震是内部位移间断边界条件，对于火山是内部压力边界条件)因此可以自然地避免内部边界的网格划分，再利用成熟的结构化网格生成技术结合网格自适应加密技术，在断层或岩浆房位置加密网格，就可以实现自动化建模和有限元计算，同时有效保证计算精度。
发明内容
[0006] 本发明正是为了解决上述技术问题而设计的一种火山和地震粘弹性形变自动建模有限元计算方法。通过等效体力替代的方式，无需划分断层/岩浆房网格，实现火山和地震引起的形变问题的自动化建模和自动化计算功能。
[0007] 本发明解决其技术问题所采用的技术方案是：
[0008] 利用成熟的结构化网格生成初始计算网格，将内边界条件用等效体力方法替换为体力项，进一步把体力项装入单元右端项，然后组装总体刚度矩阵和总体右端项，求解大型线性方程组，然后根据位移解的梯度利用网格自适应加密技术在断层/岩浆房位置加密网格，如此循环迭代2-3次，即可保证计算精度。同时采用积分型本构关系和Prony级数考虑介质的粘弹性，并且通过移动节点的方式实现地表地形起伏，用泊松方程的调和特性保证网格不产生畸变，最后利用开源的有限元库和并行求解库进行有限元并行代码编写，实现火山和地震粘弹性形变自动建模有限元计算方法。
[0009] 一种火山和地震粘弹性形变自动建模有限元计算方法，利用等效体力公式将地震剪切位错和火山膨胀压力源等效为有限元计算中的体力项，再通过局部网格自适应加密技术在地震与火山源处产生足够密的网格，省略传统地震和火山形变有限元计算中最费时的前处理步骤，自动完成建模前处理到有限元计算的全部步骤，该方法具体步骤如下：
[0010] 步骤一：地震剪切位错与火山膨胀压力源的等效体力替代
[0011] 火山膨胀压力源的等效体力替代如下：
[0012] Mogi模型代表的球状腔体压力源问题用3个正交扩张点源或3个正交张拉位错，如图1、图2和图3，这些概念上不同的源可以产生源外部相同的位移场；3对正交力偶Fh产生的位移为：
[0013] ui(x)＝-MjkGik,j(x) [1]
[0014] 其中Gik(x)是格林函数，代表源处k方向单位力产生的点接收点x处的i分量位移，并且矩张量为Mjk＝Fhδjk。代入弹性半空间内部集中力的格林函数表达式，当岩浆房半径a和岩浆房压力变化ΔP满足此时3个正交力偶可以产生等效于Mogi模型的位移场；且由3个正交力偶代表的体力项，可以写成：
[0015] ρf＝f0▽δx＝x ' [2]
[0016] 其中，δx＝x'是狄拉克函数，代表x'处(源处)的点源集中力，▽δx＝x'代表一对符号相反的脉冲，每一对力偶的强度可表示为：
[0017]
[0018] 由于狄拉克函数可用高斯分布函数的极限来近似，因此将体力项表达为高斯分布函数的形式如下：
[0019]
[0020] 其中，αxi是高斯函数在xi方向的变量，αxi取值主要取决于网格大小；式[4]中的αx1αx2αx3π3/2代表了单位体积内高斯函数的正则化系数，这种处理保证了等效体力函数足够光滑；并且，体力函数可以用来模拟不同的岩浆房形状，当αx1＝a,αx2＝b,αx3＝c且α≠b≠c时，式[5]可以用来模拟椭球状火山膨胀压力源，且椭球半轴分别位a，b和c；其他形状的压力源，可以根据弹性力学叠加原理，用式[5]的权重组合形式，即地震矩张量来表示：
[0021]
[0022] 地震剪切位错的等效体力替代如下：
[0023] 对于与地震错动等效的剪切位错，与其等效的体力项如下：
[0024]
[0025] 其中， 为等效体力，Φm表示单元K中第m个形函数，vj是位错面法向，u(r)i为位错，Cijpq为弹性系数，dΣ是单元内部位错面；
[0026] 步骤二：地震与火山源处自适应网格加密
[0027] 网格自适应加密后，会引入悬挂节点，如图4；对网格自适应加密之后悬点的处理方法如下：
[0028] 为了满足解的连续性，悬挂节点8和9的解需要满足式[7]：
[0029]
[0030] 上述u8和u9需要从有限元线性系统里消去，待求解完线性系统后通过回带获得它们的值；悬挂点的约束方程可以写成式[8]：
[0031] uiconstrained＝Cijuj [8]
[0032] 其中，
[0033] 此时，有限元线性系统，Au＝b [10]
[0034] 需要转换为：
[0035]
[0036] uiconstrained＝Cijuj [12]
[0037] 待求解式[11]后，代入式[12]得到悬挂节点的解；
[0038] 步骤三：横向非均匀椭球形地球的实现
[0039] 在计算模型中实现地形起伏和介质横向不均匀性采用不同的方法实现；对于地形起伏，首先生成无地形的球形地球，然后根据高程数据，计算出地表每个节点移动位移量，将目标节点移动到目标位置；为避免移动节点造成网格畸变，同时对地球内部节点施加移动量，内部节点移动量的确定是具体方法是将地表节点位移量作为为地球模型Poisson方程的边界条件，见式[13]，得到地球内部各节点的移动量，Poisson方程解的光滑性保证了使用这种方式添加地形不会造成太大的网格畸变；
[0040]
[0041] 对横向不均匀性以及Moho面起伏的处理，不采用移动网格的方式，而是在组建刚度矩阵时，根据Gauss积分点坐标插值得到各积分点的拉美系数和粘滞系数；
[0042] 步骤四：粘弹性本构实现
[0043] 利用玻尔兹曼叠加原理，得到maxwell体的积分型本构关系，再利用小应变假设，推导得到小应变假设下的应力递推公式：
[0044] σ(tn)＝E'ε'+σ0'
[0045] 其中，
[0046] E'＝α(Δt)πd+πV
[0047] ε'＝Δε(tn)
[0048] σ0'＝β(Δt)σd(tn-1)+σV(tn-1)。
[0049] 所述火山和地震粘弹性形变自动建模有限元计算方法，对火山膨胀压力源的体力等效，通过正交力偶和狄拉克函数的引入，得到体力项的力学表达，进一步利用高斯函数对狄拉克函数的等效，克服狄拉克函数难以在有限元方法中处理的难题，用高斯分布函数和地震矩张量的乘积得到火山膨胀压力源的体力等效项；高斯分布函数只与火山膨胀压力源的位置坐标和岩浆房三个主轴有关；地震矩张量的确定需要利用对角单位矩阵进行主轴缩放得到主轴和笛卡尔坐标系一致的椭球状火山膨胀压力源的张量表示，再进行坐标旋转得到任意主轴方向椭球状火山膨胀压力源的地震矩张量。
[0050] 所述火山和地震粘弹性形变自动建模有限元计算方法，对地震剪切位错的体力等效，先计算六面体单元被平面切割后得到的多边形，在此基础上计算地震剪切位错在多边形上的面积分，得到地震剪切位错源的等效体力项。
[0051] 所述火山和地震粘弹性形变自动建模有限元计算方法，网格自适应加密技术，采用悬挂节点的方式，这种处理方式是一种局部加密，不影响网格的整体复杂度；悬挂节点作为单元新引入的中节点，为了保证解的连续性，在有限元线性系统里面增加额外的约束方程来约束悬挂点；自适应加密进行多次迭代，但每个单元的每个面最多被细分一次，以此来保证不增加网格复杂度，即母单元细分成子单元后，新一轮的加密只能对子单元施加。
[0052] 所述火山和地震粘弹性形变自动建模有限元计算方法，基于消息传递界面的高度并行计算方法，先用主进程生成网格，然后采用METIS开源库对网格进行并行分区，对分区后的网格在各子进程上进行刚度矩阵的并行组装，然后通过进程间的通讯组装整体刚度矩阵，使用开源线性方程求解库Trillinos进行大型稀疏线性方程组的并行求解，最后各子进程将各自的计算结果进行结果文件写操作，实现后处理结果的并行显示。第六步：程序实现[0053] 由于从底层开始编写有限元程序，工作量巨大，并且调试十分繁复，因此该发明方法的有限元程序的编写主要是在C++开源有限元库deal II平台(https://www.dealii.org/)上进行。对硬件系统的要求主要是一些常用科学计算库：C++编译器、CMake编译配置工具、线性代数库BLAS和LAPACK、MPI库以支持程序并行计算、大型线性方程组求解器Trilinos和PETSc、并行网格分区库Metis、C++有限元库deal II、后处理库Paraview，这些库均是开源科学计算库，用户可以很方便的下载并安装。在编写有限元代码时，可以分为以下几个专用类，根据有限元的流程，逐一调用这些类如图5所示，就可以完成有限元代码的编写。
[0054] 利用成熟的结构化网格生成初始计算网格，将内边界条件用等效体力方法替换为体力项，进一步把体力项装入单元右端项，然后组装总体刚度矩阵和总体右端项，求解大型线性方程组，然后根据位移解的梯度利用网格自适应加密技术在断层/岩浆房位置加密网格，如此循环迭代2-3次，即可保证计算精度。同时采用积分型本构关系和Prony级数考虑介质的粘弹性，并且通过移动节点的方式实现地表地形起伏，用泊松方程的调和特性保证网格不产生畸变，最后利用开源的有限元库和并行求解库进行有限元并行代码编写，实现火山和地震粘弹性形变自动建模有限元计算方法。
[0055] 本发明的有益效果是：该发明一种火山和地震粘弹性形变自动建模有限元计算方法自动完成有限元前处理和计算，计算方法高度并行化，通过软件实现自动处理，既节省建模时间又保证计算精度。
附图说明
[0056] 图1为火山地下压力变形源的压力边界条件示意图。
[0057] 图2为火山地下压力变形源的力偶等效示意图。
[0058] 图3为火山地下压力变形源的位错等效示意图。
[0059] 图4为对网格自适应加密之后悬挂节点示意图。
[0060] 图5为有限元程序调用的各种类功能示意图。
[0061] 图6为竖直右旋断层地震剪切位错空间分布示意图。
[0062] 图7为竖直右旋断层地震剪切位错有限元与解析解对比示意图。
[0063] 图8为火山膨胀压力源有限元计算东西向位移示意图。
[0064] 图9为火山膨胀压力源有限元计算南北向位移示意图。
[0065] 图10为火山膨胀压力源有限元计算垂向位移示意图。
[0066] 图11为火山膨胀压力源有限元计算与解析解东西向位移对比示意图。
[0067] 图12为火山膨胀压力源有限元计算与解析解垂向位移对比示意图。
[0068] 图13为火山膨胀压力源有限元计算与GPS观测对比示意图。
[0069] 图14为火山膨胀压力源有限元计算与水准观测对比示意图。
具体实施方式
[0070] 下面结合附图和实施例对本发明进一步说明。
[0071] 如图1－5所示，本发明一种火山和地震粘弹性形变自动建模有限元计算方法，利用等效体力公式将地震剪切位错和火山膨胀压力源等效为有限元计算中的体力项，再通过局部网格自适应加密技术在地震与火山源处产生足够密的网格，省略传统地震和火山形变有限元计算中最费时的前处理步骤，自动完成建模前处理到有限元计算的全部步骤，该方法具体步骤如下：
[0072] 步骤一：地震剪切位错与火山膨胀压力源的等效体力替代
[0073] 火山膨胀压力源的等效体力替代如下：
[0074] Mogi模型代表的球状腔体压力源问题用3个正交扩张点源或3个正交张拉位错，如图1、图2和图3，这些概念上不同的源可以产生源外部相同的位移场；3对正交力偶Fh产生的位移为：
[0075] ui(x)＝-MjkGik,j(x) [1]
[0076] 其中Gik(x)是格林函数，代表源处k方向单位力产生的点接收点x处的i分量位移，并且矩张量为Mjk＝Fhδjk。代入弹性半空间内部集中力的格林函数表达式，当岩浆房半径a和岩浆房压力变化ΔP满足此时3个正交力偶可以产生等效于Mogi模型的位移场；且由3个正交力偶代表的体力项，可以写成：
[0077] ρf＝f0▽δx＝x ' [2]
[0078] 其中，δx＝x'是狄拉克函数，代表x'处(源处)的点源集中力，▽δx＝x'代表一对符号相反的脉冲，每一对力偶的强度可表示为：
[0079]
[0080] 由于狄拉克函数可用高斯分布函数的极限来近似，因此将体力项表达为高斯分布函数的形式如下：
[0081]
[0082] 其中，αxi是高斯函数在xi方向的变量，αxi取值主要取决于网格大小；式[4]中的αx1αx2αx3π3/2代表了单位体积内高斯函数的正则化系数，这种处理保证了等效体力函数足够光滑；并且，体力函数可以用来模拟不同的岩浆房形状，当αx1＝a,αx2＝b,αx3＝c且α≠b≠c时，式[5]可以用来模拟椭球状火山膨胀压力源，且椭球半轴分别位a，b和c；其他形状的压力源，可以根据弹性力学叠加原理，用式[5]的权重组合形式，即地震矩张量来表示：
[0083]
[0084] 地震剪切位错的等效体力替代如下：
[0085] 对于与地震错动等效的剪切位错，与其等效的体力项如下：
[0086]
[0087] 其中， 为等效体力，Φm表示单元K中第m个形函数，vj是位错面法向，u(r)i为位错，Cijpq为弹性系数，dΣ是单元内部位错面；
[0088] 步骤二：地震与火山源处自适应网格加密
[0089] 网格自适应加密后，会引入悬挂节点，如图4；对网格自适应加密之后悬点的处理方法如下：
[0090] 为了满足解的连续性，悬挂节点8和9的解需要满足式[7]：
[0091]
[0092] 上述u8和u9需要从有限元线性系统里消去，待求解完线性系统后通过回带获得它们的值；悬挂点的约束方程可以写成式[8]：
[0093] uiconstrained＝Cijuj [8]
[0094] 其中，
[0095] 此时，有限元线性系统，Au＝b [10]
[0096] 需要转换为：
[0097]
[0098] uiconstrained＝Cijuj [12]
[0099] 待求解式[11]后，代入式[12]得到悬挂节点的解；
[0100] 步骤三：横向非均匀椭球形地球的实现
[0101] 在计算模型中实现地形起伏和介质横向不均匀性采用不同的方法实现；对于地形起伏，首先生成无地形的球形地球，然后根据高程数据，计算出地表每个节点移动位移量，将目标节点移动到目标位置；为避免移动节点造成网格畸变，同时对地球内部节点施加移动量，内部节点移动量的确定是具体方法是将地表节点位移量作为为地球模型Poisson方程的边界条件，见式[13]，得到地球内部各节点的移动量，Poisson方程解的光滑性保证了使用这种方式添加地形不会造成太大的网格畸变；
[0102]
[0103] 对横向不均匀性以及Moho面起伏的处理，不采用移动网格的方式，而是在组建刚度矩阵时，根据Gauss积分点坐标插值得到各积分点的拉美系数和粘滞系数；
[0104] 步骤四：粘弹性本构实现
[0105] 利用玻尔兹曼叠加原理，得到maxwell体的积分型本构关系，再利用小应变假设，推导得到小应变假设下的应力递推公式：
[0106] σ(tn)＝E'ε'+σ0'
[0107] 其中，
[0108] E'＝α(Δt)πd+πV
[0109] ε'＝Δε(tn)
[0110] σ0'＝β(Δt)σd(tn-1)+σV(tn-1)。
[0111] 所述火山和地震粘弹性形变自动建模有限元计算方法，对火山膨胀压力源的体力等效，通过正交力偶和狄拉克函数的引入，得到体力项的力学表达，进一步利用高斯函数对狄拉克函数的等效，克服狄拉克函数难以在有限元方法中处理的难题，用高斯分布函数和地震矩张量的乘积得到火山膨胀压力源的体力等效项；高斯分布函数只与火山膨胀压力源的位置坐标和岩浆房三个主轴有关；地震矩张量的确定需要利用对角单位矩阵进行主轴缩放得到主轴和笛卡尔坐标系一致的椭球状火山膨胀压力源的张量表示，再进行坐标旋转得到任意主轴方向椭球状火山膨胀压力源的地震矩张量。
[0112] 所述火山和地震粘弹性形变自动建模有限元计算方法，对地震剪切位错的体力等效，先计算六面体单元被平面切割后得到的多边形，在此基础上计算地震剪切位错在多边形上的面积分，得到地震剪切位错源的等效体力项。
[0113] 所述火山和地震粘弹性形变自动建模有限元计算方法，网格自适应加密技术，采用悬挂节点的方式，这种处理方式是一种局部加密，不影响网格的整体复杂度；悬挂节点作为单元新引入的中节点，为了保证解的连续性，在有限元线性系统里面增加额外的约束方程来约束悬挂点；自适应加密进行多次迭代，但每个单元的每个面最多被细分一次，以此来保证不增加网格复杂度，即母单元细分成子单元后，新一轮的加密只能对子单元施加。
[0114] 所述火山和地震粘弹性形变自动建模有限元计算方法，基于消息传递界面的高度并行计算方法，先用主进程生成网格，然后采用METIS开源库对网格进行并行分区，对分区后的网格在各子进程上进行刚度矩阵的并行组装，然后通过进程间的通讯组装整体刚度矩阵，使用开源线性方程求解库Trillinos进行大型稀疏线性方程组的并行求解，最后各子进程将各自的计算结果进行结果文件写操作，实现后处理结果的并行显示。
[0115] 1.地震粘弹性形变自动化建模有限元计算应用实例
[0116] 我们计算一个竖直右旋断层如图6所示，所造成的同震变形，图7表明，位移计算结果与解析解之间很好吻合，在本文例题采用约40万单元计算时，位移结果的相对误差在大部分区域低于3％。
[0117] 2火山形变自动化建模有限元计算应用实例
[0118] 我们计算一个半径为1km，压力为10MPa，深度为5km的地下岩浆房压力变形源引起的地表形变和解析的Mogi解进行对比，见图8-图12，可以看出有限元的位移计算结果与解析解之间吻合地非常好。
[0119] 长白山火山是一座具有潜在喷发危险性的近代火山，前人对长白山火山进行了较为系统的地质、地球物理探测研究，并对长白山天池火山进行了地震活动性、形变、地球化学为主的监测研究活动，多方面研究显示长白山天池火山区下部存在低速高导岩浆房，正处于初始扰动状态向开始动荡状态的演化阶段。我们采用本发明方法，对长白山火山区进行正反演计算，含地形的网格如图10、11所示，反演得到地下岩浆房的位置和火山活动引起的地震活动之间的空间对应关系较好，得到表面速度场和现今GPS较为吻合(图13)，和水准测量的垂向位移结果也较一致(图14)。和传统的解析Mogi模型相比，计算反演的岩浆房形状为椭球状，并且充分考虑了地形，充分体现了本发明创造的优点。
[0120] 本发明不局限于上述最佳实施方式，任何人在本发明的启示下得出的其他任何与本发明相同或相近似的产品，均落在本发明的保护范围之内。
