[发明专利]非均匀椭球地球地震和地表载荷库伦应力计算方法

权利要求书

1.非均匀椭球地球地震和地表载荷库伦应力计算方法，其特征在于：利用将地震剪切位错和地表载荷分别等效为有限元计算中的体力项和面力项，再通过局部网格自适应加密技术保证计算精度，此方法可计算地震包括震间、同震和震后，以及地表载荷包括地形剥蚀、冰川消融和水库抽蓄水的加卸载引起的库伦应力变化，具体步骤如下：

步骤一：地震剪切位错与地表荷载的等效替代

包含两类情况的等效替代，对于地震剪切位错采用体力项进行等效，对于地表载荷采用面力项进行等效；

a、地震剪切位错的等效体力替代如下：

对于与地震错动等效的剪切位错，与其等效的体力项如下：

其中，为等效体力，Φ_m表示单元K中第m个形函数，v_j是位错面法向，u(r)_i为位错，C_ijpq为弹性系数，dΣ是单元内部位错面；

b、地表载荷的等效面力替代如下：

地表载荷根据获取的岩石、冰川或者水体的质量数据，和总面积，按照质量载荷除以底面积的方式，得到地表非均匀的面载荷，其中质量载荷是质量乘以重力加速度所得；

步骤二：地震剪切位错与地表荷载处网格自适应加密

网格自适应加密后，引入悬挂节点:

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

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

u_i^constrained＝C_iju_j [3]

其中，

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

需要转换为：

u_i^constrained＝C_iju_j [7]

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

步骤三：粘弹性maxwell体积分型本构递推实现

各向同性材料的三维粘弹性积分型本构关系写为：

式中，E_ijkl为松弛函数，ε_kl为应变，σ_ij为应力；在小应变假设条件下，在足够小时间增量Δt_n内应变率视为常数，经推导可得应力递推公式:

假设体积形变和剪切形变解耦，则对由1个弹簧和N个Maxwell体链并联的广义Maxwell体，体积模量和剪切模量松弛函数分别用Prony级数表示：

式中，K^∞_ijkl和G^∞_ijkl分别为时间无穷时刻的体积模量和剪切模量，N为Maxwell链数目，λ为松弛时间；认为地球介质体积模量不随时间衰减，即体积变形为弹性，剪切变形为粘弹性；对于Maxwell体，式[10]-[11]表示的剪切模量和体积模量写为：

K_ijkl(t)＝K_ijkl [12]

G_ijkl(t)＝G_ijklexp(-t/λ) [13]

体积形变和剪切形变解耦，公式[9]写为应力球张量和偏张量的形式：

σ(t_n)＝(α(Δt)π_d+π_V)Δε(t_n)+β(Δt)σ_d(t_n-1)+σ_V(t_n-1) [14]

其中，

β(Δt)＝exp(-Δt/λ)，

其中，G和K分别为剪切模量和体积模量；

将[14]改写成σ(t_n)＝E'ε'+σ₀' [16]

其中，

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

ε'＝Δε(t_n)

σ₀′＝β(Δt)σ_d(t_n-1)+σ_V(t_n-1) [17]

式[16]在形式上与带初应力的线弹性本构关系完全相同，不同的是模量和初应力随时间变化；这样，粘弹性问题就转化为带初应力的线弹性本构关系；tn时刻，用应变增量表示的虚功方程可表示为：

∫_V(δε)^TE'dV+∫_V(δε)^Tσ₀'dV-∫_V(δu)^TfdV-∫_S(δu)^TbdS＝0 [18]

通过式[18]的第3项可以按照等效体力考虑地震位错，第4项可以按照等效面力考虑地表剥蚀、冰川消融地球表面载荷过程；

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

在计算模型中添加地形起伏以及介质的横向不均匀性，二者采用不同的方法实现；

对于地形起伏，首先生成没有地形的球形地球，然后选取合适的地形数据，计算出地表每个节点需要挪动的距离，将这些节点挪动到目标位置处；为了避免挪动节点造成的网格畸变，同时对地球内部的节点施加适当的挪动量；具体方法是将地表节点的挪动量看做位移，将这个位移作为地球模型的边界条件，以地球模型为求解区域，求解一个Poisson方程，见式[19]，得到地球内部各节点需要的位移量；Poisson方程解的光滑性保证了使用这种方式添加地形不会造成太大的网格畸变；

对横向不均匀性以及Moho面起伏的处理，不采用移动网格的方式，而是在组建刚度矩阵时，根据所采用的介质参数数据，针对不同的Gauss积分点位置选取对应的参数进行计算。