[发明专利]火山和地震粘弹性形变自动建模有限元计算方法

权利要求书

1.一种火山和地震粘弹性形变自动建模有限元计算方法，其特征在于：利用等效体力公式将地震剪切位错和火山膨胀压力源等效为有限元计算中的体力项，再通过局部网格自适应加密技术在地震与火山源处产生足够密的网格，省略传统地震和火山形变有限元计算中最费时的前处理步骤，自动完成建模前处理到有限元计算的全部步骤，该方法具体步骤如下：

步骤一：地震剪切位错与火山膨胀压力源的等效体力替代

火山膨胀压力源的等效体力替代如下：

Mogi模型代表的球状腔体压力源问题用3个正交扩张点源或3个正交张拉位错，这些概念上不同的源可以产生源外部相同的位移场；3对正交力偶Fh产生的位移为：

u_i(x)＝-M_jkG_ik,j(x) [1]

其中G_ik(x)是格林函数，代表源处k方向单位力产生的点接收点x处的i分量位移，并且矩张量为M_jk＝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个形函数，v_j是位错面法向，u(r)_i为位错，C_ijpq为弹性系数，dΣ是单元内部位错面；

步骤二：地震与火山源处自适应网格加密

网格自适应加密后，会引入悬挂节点；对网格自适应加密之后悬点的处理方法如下：

为了满足解的连续性，悬挂节点8和9的解需要满足式[7]：

上述u₈和u₉需要从有限元线性系统里消去，待求解完线性系统后通过回带获得它们的值；悬挂点的约束方程可以写成式[8]：

u_i^constrained＝C_iju_j [8]

其中，

此时，有限元线性系统，Au＝b [10]

需要转换为：

u_i^constrained＝C_iju_j [12]

待求解式[11]后，代入式[12]得到悬挂节点的解；

步骤三：横向非均匀椭球形地球的实现

在计算模型中实现地形起伏和介质横向不均匀性采用不同的方法实现；对于地形起伏，首先生成无地形的球形地球，然后根据高程数据，计算出地表每个节点移动位移量，将目标节点移动到目标位置；为避免移动节点造成网格畸变，同时对地球内部节点施加移动量，内部节点移动量的确定是具体方法是将地表节点位移量作为为地球模型Poisson方程的边界条件，见式[13]，得到地球内部各节点的移动量，Poisson方程解的光滑性保证了使用这种方式添加地形不会造成太大的网格畸变；

对横向不均匀性以及Moho面起伏的处理，不采用移动网格的方式，而是在组建刚度矩阵时，根据Gauss积分点坐标插值得到各积分点的拉美系数和粘滞系数；

步骤四：粘弹性本构实现

利用玻尔兹曼叠加原理，得到maxwell体的积分型本构关系，再利用小应变假设，推导得到小应变假设下的应力递推公式：

σ(t_n)＝E'ε'+σ₀'

其中，

E'＝α(Δt)π_d+π_V

ε'＝Δε(t_n)

σ₀'＝β(Δt)σ_d(t_n-1)+σ_V(t_n-1)。