[發(fā)明專利]火山和地震粘彈性形變自動建模有限元計算方法

有效

• 申請?zhí)枺?/td>: 201811227425.0
• 申請日: 2018-10-22
• 公開（公告）號: CN109460587B
• 公開（公告）日: 2020-02-28
• 發(fā)明（設計）人: 黃祿淵;張貝;王成虎
• 申請（專利權）人: 中國地震局地殼應力研究所
• 主分類號: G06F30/23
• 分類號: G06F30/23;G06F111/10
• 代理公司: 北京興智翔達知識產權代理有限公司 11768
• 代理人: 蔣常雪

1.一種火山和地震粘彈性形變自動建模有限元計算方法，其特征在于：利用等效體力公式將地震剪切位錯和火山膨脹壓力源等效為有限元計算中的體力項，再通過局部網格自適應加密技術在地震與火山源處產生足夠密的網格，省略傳統(tǒng)地震和火山形變有限元計算中最費時的前處理步驟，自動完成建模前處理到有限元計算的全部步驟，該方法具體步驟如下：

步驟一：地震剪切位錯與火山膨脹壓力源的等效體力替代

火山膨脹壓力源的等效體力替代如下：

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Σ是單元內部位錯面；

步驟二：地震與火山源處自適應網格加密

網格自適應加密后，會引入懸掛節(jié)點；對網格自適應加密之后懸點的處理方法如下：

為了滿足解的連續(xù)性，懸掛節(jié)點8和9的解需要滿足式[7]：

上述u₈和u₉需要從有限元線性系統(tǒng)里消去，待求解完線性系統(tǒng)后通過回帶獲得它們的值；懸掛點的約束方程可以寫成式[8]：

u_i^constrained＝C_iju_j [8]

其中，

此時，有限元線性系統(tǒng)，Au＝b [10]

需要轉換為：

u_i^constrained＝C_iju_j [12]

待求解式[11]后，代入式[12]得到懸掛節(jié)點的解；

步驟三：橫向非均勻橢球形地球的實現(xiàn)

在計算模型中實現(xiàn)地形起伏和介質橫向不均勻性采用不同的方法實現(xiàn)；對于地形起伏，首先生成無地形的球形地球，然后根據高程數據，計算出地表每個節(jié)點移動位移量，將目標節(jié)點移動到目標位置；為避免移動節(jié)點造成網格畸變，同時對地球內部節(jié)點施加移動量，內部節(jié)點移動量的確定是具體方法是將地表節(jié)點位移量作為為地球模型Poisson方程的邊界條件，見式[13]，得到地球內部各節(jié)點的移動量，Poisson方程解的光滑性保證了使用這種方式添加地形不會造成太大的網格畸變；

對橫向不均勻性以及Moho面起伏的處理，不采用移動網格的方式，而是在組建剛度矩陣時，根據Gauss積分點坐標插值得到各積分點的拉美系數和粘滯系數；

步驟四：粘彈性本構實現(xiàn)

利用玻爾茲曼疊加原理，得到maxwell體的積分型本構關系，再利用小應變假設，推導得到小應變假設下的應力遞推公式：

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

其中，

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

ε'＝Δε(t_n)

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