CN109460587A - 火山和地震粘弹性形变自动建模有限元计算方法 - Google Patents
火山和地震粘弹性形变自动建模有限元计算方法 Download PDFInfo
- Publication number
- CN109460587A CN109460587A CN201811227425.0A CN201811227425A CN109460587A CN 109460587 A CN109460587 A CN 109460587A CN 201811227425 A CN201811227425 A CN 201811227425A CN 109460587 A CN109460587 A CN 109460587A
- Authority
- CN
- China
- Prior art keywords
- finite element
- volcano
- seismic
- grid
- earthquake
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 46
- CCEKAJIANROZEO-UHFFFAOYSA-N sulfluramid Chemical group CCNS(=O)(=O)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)C(F)(F)F CCEKAJIANROZEO-UHFFFAOYSA-N 0.000 title claims abstract description 42
- 238000004364 calculation method Methods 0.000 claims abstract description 74
- 238000005516 engineering process Methods 0.000 claims abstract description 11
- 238000010008 shearing Methods 0.000 claims abstract description 10
- 230000003044 adaptive effect Effects 0.000 claims abstract description 9
- 238000006073 displacement reaction Methods 0.000 claims description 31
- 239000000725 suspension Substances 0.000 claims description 19
- 239000011159 matrix material Substances 0.000 claims description 13
- 238000005315 distribution function Methods 0.000 claims description 12
- 230000008569 process Effects 0.000 claims description 12
- 238000012545 processing Methods 0.000 claims description 10
- 238000012876 topography Methods 0.000 claims description 10
- 230000014509 gene expression Effects 0.000 claims description 8
- 230000010354 integration Effects 0.000 claims description 6
- 239000002002 slurry Substances 0.000 claims description 6
- 241001191378 Moho Species 0.000 claims description 5
- 239000011435 rock Substances 0.000 claims description 5
- 230000008859 change Effects 0.000 claims description 4
- 238000005192 partition Methods 0.000 claims description 4
- 238000012805 post-processing Methods 0.000 claims description 4
- 238000006467 substitution reaction Methods 0.000 claims description 4
- 241001229889 Metis Species 0.000 claims description 3
- 230000005540 biological transmission Effects 0.000 claims description 3
- 238000004891 communication Methods 0.000 claims description 3
- 238000004613 tight binding model Methods 0.000 claims description 3
- 238000002203 pretreatment Methods 0.000 abstract description 4
- 230000009286 beneficial effect Effects 0.000 abstract description 2
- 238000007781 pre-processing Methods 0.000 abstract description 2
- 238000003325 tomography Methods 0.000 abstract 1
- 230000006870 function Effects 0.000 description 22
- 238000010586 diagram Methods 0.000 description 14
- 230000000694 effects Effects 0.000 description 5
- 238000011160 research Methods 0.000 description 5
- 238000005259 measurement Methods 0.000 description 2
- 238000000053 physical method Methods 0.000 description 2
- 230000008961 swelling Effects 0.000 description 2
- 235000010724 Wisteria floribunda Nutrition 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000019771 cognition Effects 0.000 description 1
- 230000001066 destructive effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 238000012502 risk assessment Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 230000004083 survival effect Effects 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种火山和地震粘弹性形变自动建模有限元计算方法，利用等效体力公式将地震剪切位错和火山膨胀压力源等效为有限元计算中的体力项，再通过局部网格自适应加密技术在地震与火山源处产生足够密的网格，即可省略传统地震和火山形变有限元计算中最费时的前处理步骤‑‑划分断层/火山岩浆房网格，实现火山和地震粘弹性形变的自动建模功能。本发明的有益效果是：该发明一种火山和地震粘弹性形变自动建模有限元计算方法自动完成有限元前处理和计算，计算方法高度并行化，通过软件实现自动处理，既节省建模时间又保证计算精度。
Description
技术领域
本发明涉及地球物理领域，具体地说，涉及一种火山和地震粘弹性形变自动建模有限元计算方法。主要应用于测震学和火山形变学领域，其功能是在无需划出断层和岩浆房网格的情况下，利用有限元计算因为地震静态滑动分布/火山膨胀压力源引起的形变以及应力变化，完全实现自动化，节省大量前处理建模时间。
背景技术
地震、火山是地球科学的重要研究对象，二者带来的巨大破坏力都对人类生存具有潜在威胁，但目前相关的解析计算方法仍存在较大局限性：大地震往往发生于俯冲带地区或者板内的地形和介质突变区，比如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面起伏、介质非均匀性，在地震和火山应变、应力计算中尽可能模拟真实地球，对地震危险性分析和火山区岩浆活动的认知具有重要意义。
但目前通用的地震和火山形变的有限元计算方法仍然存在改进空间，最突出的问题就是前处理网格划分十分繁琐。在传统有限元方法里，地震破裂的模拟往往需要在前处理中划出具体的断层破裂面，而火山地区的地下压力膨胀源一般要在前处理中刻画出球/椭球的岩浆房网格然后施加面力边界条件。但这种前处理方式存在以下困难：1、当断层数量较少、且断层主要为高倾角或者近直立时，前处理建模难度相对较低，但当断层倾角较低或者断层数量较多时，前处理往往十分费时且网格质量不能保证；2、当岩浆房形状复杂、或者岩浆房数量较多时，前处理建模难度大，且往往只能采用四面体网格，计算精度低于六面体网格。
鉴于以上实际应用原因，非常有必要发展一种新的火山和地震粘弹性形变有限元计算方法，要求这种计算方法自动完成有限元前处理和计算，既节省建模时间又保证计算精度。
如果还按照传统有限元方法，寻找一种策略能让程序自动划分出断层和岩浆房网格，显然这种策略的难度很高，而且难以保证网格不过度畸形。因此，我们采用另一种策略，即把断层(剪切位错)和岩浆房膨胀(膨胀位错)等效为体力。此方法不需要像传统有限元设置内边界条件(对于地震是内部位移间断边界条件，对于火山是内部压力边界条件)因此可以自然地避免内部边界的网格划分，再利用成熟的结构化网格生成技术结合网格自适应加密技术，在断层或岩浆房位置加密网格，就可以实现自动化建模和有限元计算，同时有效保证计算精度。
发明内容
本发明正是为了解决上述技术问题而设计的一种火山和地震粘弹性形变自动建模有限元计算方法。通过等效体力替代的方式，无需划分断层/岩浆房网格，实现火山和地震引起的形变问题的自动化建模和自动化计算功能。
本发明解决其技术问题所采用的技术方案是：
利用成熟的结构化网格生成初始计算网格，将内边界条件用等效体力方法替换为体力项，进一步把体力项装入单元右端项，然后组装总体刚度矩阵和总体右端项，求解大型线性方程组，然后根据位移解的梯度利用网格自适应加密技术在断层/岩浆房位置加密网格，如此循环迭代2-3次，即可保证计算精度。同时采用积分型本构关系和Prony级数考虑介质的粘弹性，并且通过移动节点的方式实现地表地形起伏，用泊松方程的调和特性保证网格不产生畸变，最后利用开源的有限元库和并行求解库进行有限元并行代码编写，实现火山和地震粘弹性形变自动建模有限元计算方法。
一种火山和地震粘弹性形变自动建模有限元计算方法，利用等效体力公式将地震剪切位错和火山膨胀压力源等效为有限元计算中的体力项，再通过局部网格自适应加密技术在地震与火山源处产生足够密的网格，省略传统地震和火山形变有限元计算中最费时的前处理步骤，自动完成建模前处理到有限元计算的全部步骤，该方法具体步骤如下：
步骤一：地震剪切位错与火山膨胀压力源的等效体力替代
火山膨胀压力源的等效体力替代如下：
Mogi模型代表的球状腔体压力源问题用3个正交扩张点源或3个正交张拉位错，如图1、图2和图3，这些概念上不同的源可以产生源外部相同的位移场； 3对正交力偶Fh产生的位移为：
ui(x)＝-MjkGik,j(x) [1]
其中Gik(x)是格林函数，代表源处k方向单位力产生的点接收点x处的i分量位移，并且矩张量为Mjk＝Fhδjk。代入弹性半空间内部集中力的格林函数表达式，当岩浆房半径a和岩浆房压力变化ΔP满足此时3个正交力偶可以产生等效于Mogi模型的位移场；且由3个正交力偶代表的体力项，可以写成：
ρf＝f0▽δx＝x' [2]
其中，δ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Σ是单元内部位错面；
步骤二：地震与火山源处自适应网格加密
网格自适应加密后，会引入悬挂节点，如图4；对网格自适应加密之后悬点的处理方法如下：
为了满足解的连续性，悬挂节点8和9的解需要满足式[7]：
上述u8和u9需要从有限元线性系统里消去，待求解完线性系统后通过回带获得它们的值；悬挂点的约束方程可以写成式[8]：
ui constrained＝Cijuj [8]
其中，
此时，有限元线性系统，Au＝b [10]
需要转换为：
ui constrained＝Cijuj [12]
待求解式[11]后，代入式[12]得到悬挂节点的解；
步骤三：横向非均匀椭球形地球的实现
在计算模型中实现地形起伏和介质横向不均匀性采用不同的方法实现；对于地形起伏，首先生成无地形的球形地球，然后根据高程数据，计算出地表每个节点移动位移量，将目标节点移动到目标位置；为避免移动节点造成网格畸变，同时对地球内部节点施加移动量，内部节点移动量的确定是具体方法是将地表节点位移量作为为地球模型Poisson方程的边界条件，见式[13]，得到地球内部各节点的移动量，Poisson方程解的光滑性保证了使用这种方式添加地形不会造成太大的网格畸变；
对横向不均匀性以及Moho面起伏的处理，不采用移动网格的方式，而是在组建刚度矩阵时，根据Gauss积分点坐标插值得到各积分点的拉美系数和粘滞系数；
步骤四：粘弹性本构实现
利用玻尔兹曼叠加原理，得到maxwell体的积分型本构关系，再利用小应变假设，推导得到小应变假设下的应力递推公式：
σ(tn)＝E'ε'+σ0'
其中，
E'＝α(Δt)πd+πV
ε'＝Δε(tn)
σ0'＝β(Δt)σd(tn-1)+σV(tn-1)。
所述火山和地震粘弹性形变自动建模有限元计算方法，对火山膨胀压力源的体力等效，通过正交力偶和狄拉克函数的引入，得到体力项的力学表达，进一步利用高斯函数对狄拉克函数的等效，克服狄拉克函数难以在有限元方法中处理的难题，用高斯分布函数和地震矩张量的乘积得到火山膨胀压力源的体力等效项；高斯分布函数只与火山膨胀压力源的位置坐标和岩浆房三个主轴有关；地震矩张量的确定需要利用对角单位矩阵进行主轴缩放得到主轴和笛卡尔坐标系一致的椭球状火山膨胀压力源的张量表示，再进行坐标旋转得到任意主轴方向椭球状火山膨胀压力源的地震矩张量。
所述火山和地震粘弹性形变自动建模有限元计算方法，对地震剪切位错的体力等效，先计算六面体单元被平面切割后得到的多边形，在此基础上计算地震剪切位错在多边形上的面积分，得到地震剪切位错源的等效体力项。
所述火山和地震粘弹性形变自动建模有限元计算方法，网格自适应加密技术，采用悬挂节点的方式，这种处理方式是一种局部加密，不影响网格的整体复杂度；悬挂节点作为单元新引入的中节点，为了保证解的连续性，在有限元线性系统里面增加额外的约束方程来约束悬挂点；自适应加密进行多次迭代，但每个单元的每个面最多被细分一次，以此来保证不增加网格复杂度，即母单元细分成子单元后，新一轮的加密只能对子单元施加。
所述火山和地震粘弹性形变自动建模有限元计算方法，基于消息传递界面的高度并行计算方法，先用主进程生成网格，然后采用METIS开源库对网格进行并行分区，对分区后的网格在各子进程上进行刚度矩阵的并行组装，然后通过进程间的通讯组装整体刚度矩阵，使用开源线性方程求解库Trillinos进行大型稀疏线性方程组的并行求解，最后各子进程将各自的计算结果进行结果文件写操作，实现后处理结果的并行显示。第六步：程序实现
由于从底层开始编写有限元程序，工作量巨大，并且调试十分繁复，因此该发明方法的有限元程序的编写主要是在C++开源有限元库deal II平台 (https://www.dealii.org/)上进行。对硬件系统的要求主要是一些常用科学计算库：C++编译器、CMake编译配置工具、线性代数库BLAS和LAPACK、MPI 库以支持程序并行计算、大型线性方程组求解器Trilinos和PETSc、并行网格分区库Metis、C++有限元库deal II、后处理库Paraview，这些库均是开源科学计算库，用户可以很方便的下载并安装。在编写有限元代码时，可以分为以下几个专用类，根据有限元的流程，逐一调用这些类如图5所示，就可以完成有限元代码的编写。
利用成熟的结构化网格生成初始计算网格，将内边界条件用等效体力方法替换为体力项，进一步把体力项装入单元右端项，然后组装总体刚度矩阵和总体右端项，求解大型线性方程组，然后根据位移解的梯度利用网格自适应加密技术在断层/岩浆房位置加密网格，如此循环迭代2-3次，即可保证计算精度。同时采用积分型本构关系和Prony级数考虑介质的粘弹性，并且通过移动节点的方式实现地表地形起伏，用泊松方程的调和特性保证网格不产生畸变，最后利用开源的有限元库和并行求解库进行有限元并行代码编写，实现火山和地震粘弹性形变自动建模有限元计算方法。
本发明的有益效果是：该发明一种火山和地震粘弹性形变自动建模有限元计算方法自动完成有限元前处理和计算，计算方法高度并行化，通过软件实现自动处理，既节省建模时间又保证计算精度。
附图说明
图1为火山地下压力变形源的压力边界条件示意图。
图2为火山地下压力变形源的力偶等效示意图。
图3为火山地下压力变形源的位错等效示意图。
图4为对网格自适应加密之后悬挂节点示意图。
图5为有限元程序调用的各种类功能示意图。
图6为竖直右旋断层地震剪切位错空间分布示意图。
图7为竖直右旋断层地震剪切位错有限元与解析解对比示意图。
图8为火山膨胀压力源有限元计算东西向位移示意图。
图9为火山膨胀压力源有限元计算南北向位移示意图。
图10为火山膨胀压力源有限元计算垂向位移示意图。
图11为火山膨胀压力源有限元计算与解析解东西向位移对比示意图。
图12为火山膨胀压力源有限元计算与解析解垂向位移对比示意图。
图13为火山膨胀压力源有限元计算与GPS观测对比示意图。
图14为火山膨胀压力源有限元计算与水准观测对比示意图。
具体实施方式
下面结合附图和实施例对本发明进一步说明。
如图1－5所示，本发明一种火山和地震粘弹性形变自动建模有限元计算方法，利用等效体力公式将地震剪切位错和火山膨胀压力源等效为有限元计算中的体力项，再通过局部网格自适应加密技术在地震与火山源处产生足够密的网格，省略传统地震和火山形变有限元计算中最费时的前处理步骤，自动完成建模前处理到有限元计算的全部步骤，该方法具体步骤如下：
步骤一：地震剪切位错与火山膨胀压力源的等效体力替代
火山膨胀压力源的等效体力替代如下：
Mogi模型代表的球状腔体压力源问题用3个正交扩张点源或3个正交张拉位错，如图1、图2和图3，这些概念上不同的源可以产生源外部相同的位移场； 3对正交力偶Fh产生的位移为：
ui(x)＝-MjkGik,j(x) [1]
其中Gik(x)是格林函数，代表源处k方向单位力产生的点接收点x处的i分量位移，并且矩张量为Mjk＝Fhδjk。代入弹性半空间内部集中力的格林函数表达式，当岩浆房半径a和岩浆房压力变化ΔP满足此时3个正交力偶可以产生等效于Mogi模型的位移场；且由3个正交力偶代表的体力项，可以写成：
ρf＝f0▽δx＝x' [2]
其中，δ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Σ是单元内部位错面；
步骤二：地震与火山源处自适应网格加密
网格自适应加密后，会引入悬挂节点，如图4；对网格自适应加密之后悬点的处理方法如下：
为了满足解的连续性，悬挂节点8和9的解需要满足式[7]：
上述u8和u9需要从有限元线性系统里消去，待求解完线性系统后通过回带获得它们的值；悬挂点的约束方程可以写成式[8]：
ui constrained＝Cijuj [8]
其中，
此时，有限元线性系统，Au＝b [10]
需要转换为：
ui constrained＝Cijuj [12]
待求解式[11]后，代入式[12]得到悬挂节点的解；
步骤三：横向非均匀椭球形地球的实现
在计算模型中实现地形起伏和介质横向不均匀性采用不同的方法实现；对于地形起伏，首先生成无地形的球形地球，然后根据高程数据，计算出地表每个节点移动位移量，将目标节点移动到目标位置；为避免移动节点造成网格畸变，同时对地球内部节点施加移动量，内部节点移动量的确定是具体方法是将地表节点位移量作为为地球模型Poisson方程的边界条件，见式[13]，得到地球内部各节点的移动量，Poisson方程解的光滑性保证了使用这种方式添加地形不会造成太大的网格畸变；
对横向不均匀性以及Moho面起伏的处理，不采用移动网格的方式，而是在组建刚度矩阵时，根据Gauss积分点坐标插值得到各积分点的拉美系数和粘滞系数；
步骤四：粘弹性本构实现
利用玻尔兹曼叠加原理，得到maxwell体的积分型本构关系，再利用小应变假设，推导得到小应变假设下的应力递推公式：
σ(tn)＝E'ε'+σ0'
其中，
E'＝α(Δt)πd+πV
ε'＝Δε(tn)
σ0'＝β(Δt)σd(tn-1)+σV(tn-1)。
所述火山和地震粘弹性形变自动建模有限元计算方法，对火山膨胀压力源的体力等效，通过正交力偶和狄拉克函数的引入，得到体力项的力学表达，进一步利用高斯函数对狄拉克函数的等效，克服狄拉克函数难以在有限元方法中处理的难题，用高斯分布函数和地震矩张量的乘积得到火山膨胀压力源的体力等效项；高斯分布函数只与火山膨胀压力源的位置坐标和岩浆房三个主轴有关；地震矩张量的确定需要利用对角单位矩阵进行主轴缩放得到主轴和笛卡尔坐标系一致的椭球状火山膨胀压力源的张量表示，再进行坐标旋转得到任意主轴方向椭球状火山膨胀压力源的地震矩张量。
所述火山和地震粘弹性形变自动建模有限元计算方法，对地震剪切位错的体力等效，先计算六面体单元被平面切割后得到的多边形，在此基础上计算地震剪切位错在多边形上的面积分，得到地震剪切位错源的等效体力项。
所述火山和地震粘弹性形变自动建模有限元计算方法，网格自适应加密技术，采用悬挂节点的方式，这种处理方式是一种局部加密，不影响网格的整体复杂度；悬挂节点作为单元新引入的中节点，为了保证解的连续性，在有限元线性系统里面增加额外的约束方程来约束悬挂点；自适应加密进行多次迭代，但每个单元的每个面最多被细分一次，以此来保证不增加网格复杂度，即母单元细分成子单元后，新一轮的加密只能对子单元施加。
所述火山和地震粘弹性形变自动建模有限元计算方法，基于消息传递界面的高度并行计算方法，先用主进程生成网格，然后采用METIS开源库对网格进行并行分区，对分区后的网格在各子进程上进行刚度矩阵的并行组装，然后通过进程间的通讯组装整体刚度矩阵，使用开源线性方程求解库Trillinos进行大型稀疏线性方程组的并行求解，最后各子进程将各自的计算结果进行结果文件写操作，实现后处理结果的并行显示。
1.地震粘弹性形变自动化建模有限元计算应用实例
我们计算一个竖直右旋断层如图6所示，所造成的同震变形，图7表明，位移计算结果与解析解之间很好吻合，在本文例题采用约40万单元计算时，位移结果的相对误差在大部分区域低于3％。
2火山形变自动化建模有限元计算应用实例
我们计算一个半径为1km，压力为10MPa，深度为5km的地下岩浆房压力变形源引起的地表形变和解析的Mogi解进行对比，见图8-图12，可以看出有限元的位移计算结果与解析解之间吻合地非常好。
长白山火山是一座具有潜在喷发危险性的近代火山，前人对长白山火山进行了较为系统的地质、地球物理探测研究，并对长白山天池火山进行了地震活动性、形变、地球化学为主的监测研究活动，多方面研究显示长白山天池火山区下部存在低速高导岩浆房，正处于初始扰动状态向开始动荡状态的演化阶段。我们采用本发明方法，对长白山火山区进行正反演计算，含地形的网格如图10、 11所示，反演得到地下岩浆房的位置和火山活动引起的地震活动之间的空间对应关系较好，得到表面速度场和现今GPS较为吻合(图13)，和水准测量的垂向位移结果也较一致(图14)。和传统的解析Mogi模型相比，计算反演的岩浆房形状为椭球状，并且充分考虑了地形，充分体现了本发明创造的优点。
本发明不局限于上述最佳实施方式，任何人在本发明的启示下得出的其他任何与本发明相同或相近似的产品，均落在本发明的保护范围之内。
Claims (5)
1.一种火山和地震粘弹性形变自动建模有限元计算方法，其特征在于：利用等效体力公式将地震剪切位错和火山膨胀压力源等效为有限元计算中的体力项，再通过局部网格自适应加密技术在地震与火山源处产生足够密的网格，省略传统地震和火山形变有限元计算中最费时的前处理步骤，自动完成建模前处理到有限元计算的全部步骤，该方法具体步骤如下：
步骤一：地震剪切位错与火山膨胀压力源的等效体力替代
火山膨胀压力源的等效体力替代如下：
Mogi模型代表的球状腔体压力源问题用3个正交扩张点源或3个正交张拉位错，这些概念上不同的源可以产生源外部相同的位移场；3对正交力偶Fh产生的位移为：
ui(x)＝-MjkGik,j(x) [1]
其中Gik(x)是格林函数，代表源处k方向单位力产生的点接收点x处的i分量位移，并且矩张量为Mjk＝Fhδjk；代入弹性半空间内部集中力的格林函数表达式，当岩浆房半径a和岩浆房压力变化ΔP满足此时3个正交力偶可以产生等效于Mogi模型的位移场；且由3个正交力偶代表的体力项，可以写成：
其中，δ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]：
ui constrained＝Cijuj [8]
其中，
此时，有限元线性系统，Au＝b [10]
需要转换为：
ui constrained＝Cijuj [12]
待求解式[11]后，代入式[12]得到悬挂节点的解；
步骤三：横向非均匀椭球形地球的实现
在计算模型中实现地形起伏和介质横向不均匀性采用不同的方法实现；对于地形起伏，首先生成无地形的球形地球，然后根据高程数据，计算出地表每个节点移动位移量，将目标节点移动到目标位置；为避免移动节点造成网格畸变，同时对地球内部节点施加移动量，内部节点移动量的确定是具体方法是将地表节点位移量作为为地球模型Poisson方程的边界条件，见式[13]，得到地球内部各节点的移动量，Poisson方程解的光滑性保证了使用这种方式添加地形不会造成太大的网格畸变；
对横向不均匀性以及Moho面起伏的处理，不采用移动网格的方式，而是在组建刚度矩阵时，根据Gauss积分点坐标插值得到各积分点的拉美系数和粘滞系数；
步骤四：粘弹性本构实现
利用玻尔兹曼叠加原理，得到maxwell体的积分型本构关系，再利用小应变假设，推导得到小应变假设下的应力递推公式：
σ(tn)＝E'ε'+σ0'
其中，
E'＝α(Δt)πd+πV
ε'＝Δε(tn)
σ0'＝β(Δt)σd(tn-1)+σV(tn-1)。
2.根据权利要求1所述的火山和地震粘弹性形变自动建模有限元计算方法，其特征在于：对火山膨胀压力源的体力等效，通过正交力偶和狄拉克函数的引入，得到体力项的力学表达，进一步利用高斯函数对狄拉克函数的等效，克服狄拉克函数难以在有限元方法中处理的难题，用高斯分布函数和地震矩张量的乘积得到火山膨胀压力源的体力等效项；高斯分布函数只与火山膨胀压力源的位置坐标和岩浆房三个主轴有关；地震矩张量的确定需要利用对角单位矩阵进行主轴缩放得到主轴和笛卡尔坐标系一致的椭球状火山膨胀压力源的张量表示，再进行坐标旋转得到任意主轴方向椭球状火山膨胀压力源的地震矩张量。
3.根据权利要求1所述的火山和地震粘弹性形变自动建模有限元计算方法，其特征在于：对地震剪切位错的体力等效，先计算六面体单元被平面切割后得到的多边形，在此基础上计算地震剪切位错在多边形上的面积分，得到地震剪切位错源的等效体力项。
4.根据权利要求1所述的火山和地震粘弹性形变自动建模有限元计算方法，其特征在于：网格自适应加密技术，采用悬挂节点的方式，这种处理方式是一种局部加密，不影响网格的整体复杂度；悬挂节点作为单元新引入的中节点，为了保证解的连续性，在有限元线性系统里面增加额外的约束方程来约束悬挂点；自适应加密进行多次迭代，但每个单元的每个面最多被细分一次，以此来保证不增加网格复杂度，即母单元细分成子单元后，新一轮的加密只能对子单元施加。
5.根据权利要求1所述的火山和地震粘弹性形变自动建模有限元计算方法，其特征在于：基于消息传递界面的高度并行计算方法，先用主进程生成网格，然后采用METIS开源库对网格进行并行分区，对分区后的网格在各子进程上进行刚度矩阵的并行组装，然后通过进程间的通讯组装整体刚度矩阵，使用开源线性方程求解库Trillinos进行大型稀疏线性方程组的并行求解，最后各子进程将各自的计算结果进行结果文件写操作，实现后处理结果的并行显示。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811227425.0A CN109460587B (zh) | 2018-10-22 | 2018-10-22 | 火山和地震粘弹性形变自动建模有限元计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811227425.0A CN109460587B (zh) | 2018-10-22 | 2018-10-22 | 火山和地震粘弹性形变自动建模有限元计算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109460587A true CN109460587A (zh) | 2019-03-12 |
| CN109460587B CN109460587B (zh) | 2020-02-28 |
Family
ID=65608026
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811227425.0A Active CN109460587B (zh) | 2018-10-22 | 2018-10-22 | 火山和地震粘弹性形变自动建模有限元计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109460587B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113031109A (zh) * | 2019-12-25 | 2021-06-25 | 中国石油天然气股份有限公司 | 岩石物理建模方法及装置 |
| CN117934747A (zh) * | 2024-03-22 | 2024-04-26 | 山东省地震工程研究院 | 一种基于激光点云数据的活动构造地貌三维模型构建方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20100042379A1 (en) * | 2007-03-20 | 2010-02-18 | Karel Minnaar | Framework to Determine the Capacity of A Structure |
| CN103279991A (zh) * | 2013-04-16 | 2013-09-04 | 西南石油大学 | 一种利用离散裂缝端点变形数模提高油藏开发效果的方法 |
| CN103913772A (zh) * | 2014-04-02 | 2014-07-09 | 西南石油大学 | 基于储层地质力学参数的微地震事件正演模拟方法 |
| CN106646597A (zh) * | 2016-12-14 | 2017-05-10 | 中国石油大学(北京) | 基于弹簧网络模型的正演模拟方法及装置 |
| CN107798156A (zh) * | 2016-09-02 | 2018-03-13 | 赵建国 | 一种频率域2.5维粘弹性波数值模拟方法及装置 |
-
2018
- 2018-10-22 CN CN201811227425.0A patent/CN109460587B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20100042379A1 (en) * | 2007-03-20 | 2010-02-18 | Karel Minnaar | Framework to Determine the Capacity of A Structure |
| CN103279991A (zh) * | 2013-04-16 | 2013-09-04 | 西南石油大学 | 一种利用离散裂缝端点变形数模提高油藏开发效果的方法 |
| CN103913772A (zh) * | 2014-04-02 | 2014-07-09 | 西南石油大学 | 基于储层地质力学参数的微地震事件正演模拟方法 |
| CN107798156A (zh) * | 2016-09-02 | 2018-03-13 | 赵建国 | 一种频率域2.5维粘弹性波数值模拟方法及装置 |
| CN106646597A (zh) * | 2016-12-14 | 2017-05-10 | 中国石油大学(北京) | 基于弹簧网络模型的正演模拟方法及装置 |
Non-Patent Citations (2)
| Title |
|---|
| 张贝等: "有限元模拟弹性位错的等效体力方法", 《地球物理学报》 * |
| 黄禄渊等: "地震位错模拟在有限元软件MSC.Marc中的实现", 《华北地震科学》 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113031109A (zh) * | 2019-12-25 | 2021-06-25 | 中国石油天然气股份有限公司 | 岩石物理建模方法及装置 |
| CN113031109B (zh) * | 2019-12-25 | 2024-02-02 | 中国石油天然气股份有限公司 | 岩石物理建模方法及装置 |
| CN117934747A (zh) * | 2024-03-22 | 2024-04-26 | 山东省地震工程研究院 | 一种基于激光点云数据的活动构造地貌三维模型构建方法 |
| CN117934747B (zh) * | 2024-03-22 | 2024-06-07 | 山东省地震工程研究院 | 一种基于激光点云数据的活动构造地貌三维模型构建方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109460587B (zh) | 2020-02-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Dalguer et al. | Simulation of tensile crack generation by three‐dimensional dynamic shear rupture propagation during an earthquake | |
| Hasan et al. | Minimum fault coverage in reconfigurable arrays | |
| Rietmann et al. | Forward and adjoint simulations of seismic wave propagation on emerging large-scale GPU architectures | |
| CA2816815C (en) | Systems and methods for generating updates of geological models | |
| Borker et al. | Mesh adaptation framework for embedded boundary methods for computational fluid dynamics and fluid‐structure interaction | |
| CN105467443A (zh) | 一种三维各向异性弹性波数值模拟方法及系统 | |
| Maerten | Adaptive cross-approximation applied to the solution of system of equations and post-processing for 3D elastostatic problems using the boundary element method | |
| CN109460587B (zh) | 火山和地震粘弹性形变自动建模有限元计算方法 | |
| Beyabanaki et al. | Applying disk-based discontinuous deformation analysis (DDA) to simulate Donghekou landslide triggered by the Wenchuan earthquake | |
| Xue et al. | An efficient GPU implementation for locating micro-seismic sources using 3D elastic wave time-reversal imaging | |
| CN114139335A (zh) | 基于单松弛时间格子玻尔兹曼模型的粘滞声波模拟方法 | |
| Muratov et al. | Grid-characteristic method on unstructured tetrahedral meshes | |
| CN109490948B (zh) | 地震声学波动方程矢量并行计算方法 | |
| CN106662665B (zh) | 用于更快速的交错网格处理的重新排序的插值和卷积 | |
| Nurtas et al. | Convolutional neural networks as a method to solve estimation problem of acoustic wave propagation in poroelastic media | |
| Nurtas et al. | 2-D Finite Element method using" eScript" for acoustic wave propagation | |
| Brunsvik et al. | Three-dimensional paganica fault morphology obtained from hypocenter clustering (L'Aquila 2009 seismic sequence, Central Italy) | |
| CN109270590B (zh) | 非均匀椭球地球地震和地表载荷库伦应力计算方法 | |
| Shvarev et al. | The method to model microseismic events during hydrofracture propagation | |
| Landinez et al. | First steps on modelling wave propagation in isotropic-heterogeneous media: Numerical simulation of P–SV waves | |
| Morra et al. | Fresh Outlook in Numerical Methods for Geodynamics | |
| Singh et al. | GPU-based 3D anisotropic elastic modeling using mimetic finite differences | |
| Assonitis et al. | A shock-fitting technique for 2D/3D flows with interactions using structured grids | |
| Breil et al. | Second-order extension in space and time for a 3D cell-centered Lagrangian scheme | |
| Assonitis et al. | Numerical simulations of shock interactions on 3D structured grids using a shock-fitting approach |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CP03 | Change of name, title or address |
Address after: 100085, Anning Road, Beijing, Haidian District, No. 1 Patentee after: National natural disaster prevention and Control Research Institute, Ministry of emergency management Address before: 100085 No. 1 Anning Road, Xisanqi, Haidian District, Beijing. Patentee before: THE INSTITUTE OF CRUSTAL DYNAMICS, CHINA EARTHQUAKE ADMINISTRATION |
|
| CP03 | Change of name, title or address |