無限元-譜元混合法在2.5維引力位計算中的應用*

任駿聲,張懷

(中國科學院大學地球與行星科學學院 中國科學院計算地球動力學重點實驗室, 北京 100049) (2022年3月14日收稿; 2022年4月15日收修改稿)

引力位從地球內部到地球外部無處不在。當大地震發生時,地球內部物質產生位移,引起引力位擾動。因此,在地球自由振蕩以及低頻大尺度(如全球尺度)地震波場數值模擬中,引力位計算是一個不可避免的問題。但無論引力位還是引力位擾動,在地球內部均滿足泊松方程,在地球外滿足拉普拉斯方程,其邊界條件均為在無窮遠處衰減為零[1],這迫使我們需要利用數值方法求解一個無界微分方程。長期以來,為了避免求解該無界微分方程而采用Cowling近似[2],舍棄掉引力位擾動項。這會造成模擬波形的不準確,從而導致精細頻譜分析困難,特別是對于地球自由振蕩[3]、固體潮汐[4]等問題。

針對此類問題,常規的處理方法是將求解域設定在有限的范圍內。例如將地球內部網格向外延伸到足夠遠的地球外(如數倍地球半徑處),然后在這一人為指定的邊界處設置邊界條件,只要延伸得足夠遠,則對地球內部及近地表的引力位影響較小。該邊界條件的選取有多種方法：例如Cai和Wang[5]提出的r-4近似假設,用于近似表示位置r處的函數值;Chaljub和Valette[6]采用球諧函數展開構建DtN算子,用于該邊界條件的處理;Xu和Huo[7]使用Legendre多項式近似表示邊界處的值。這些處理總體上都是通過對邊界條件取值的近似來滿足所需的計算精度。事實上,對于無界泊松問題,數學上很早就有研究[8],只是一些方法沒有在地學中得到使用。

本文在延拓網格的基礎上,運用無限元映射[9]的思想,直接通過構造插值多項式,自然地滿足無窮遠趨于零的邊界條件,從而避免了對邊界條件取值的近似處理,因而結果更準確。

本文采用譜元法[10]作為數值求解方法。譜元法是一種綜合型方法,既保留了傳統有限元方法處理復雜網格和邊界的便利性,又確保了類似譜配置方法的高精度?？紤]到全球尺度模擬的計算量和模型的復雜度(如S40RTS[11]地幔模型),本文采用2.5維[12-13]的控制方程。在2.5維假設下,引力位的計算簡化到二維D形區域。由于對稱軸的原因,軸單元的數值積分需要進行改進。本文采用Gerritsma和Phillips[14]提出的歸一化方法,在軸單元引入GRL(Gauss-Radau-Legendre)數值積分。

為分析此模擬方法的有效性,本文開展了數值試驗,以三維球對稱實心球為例,數值模擬其引力位,并與解析解做對比。

1 方法

1.1 控制方程

對于引力位Ψ以及引力位擾動φ,均滿足泊松-拉普拉斯方程[1]：

(1)

(2)

其中：g為地球的引力常數,取(6.673 84±0.000 12)×10-11m3·kg-1·s-2;ρ為地球密度;d為地震引起的位移變化;V為整個地球。

為簡化問題,方程(1)、(2)可以簡寫為

(3)

其中s代表源項,是位置x的函數,在地球外部取0,在地球內部與密度有關。方程(3)對應的邊界條件為

(4)

采用有限元方法,對方程(3)做弱形式有

(5)

其中ψ為測試函數。

1.2 離散與網格剖分

在柱坐標系(r,θ,z)下,對3維方程(5)進行簡化(2.5維)?？紤]對稱性,關于θ的運算可以忽略,此時位置x(r,θ,z)變為x(r,z)。相應地,體積元dV=rdrdθdz變成

dΩ=rdrdz.

(6)

按照譜元法的要求,對求解域Ω進行四邊形網格剖分,如圖1藍色部分所示。本文采用doubling算法,在介質間斷面處做加密處理。對于地球外部的網格延拓方案：首先,以地球網格為基礎,對近地表處的網格進行局部加密;然后沿r向無窮遠方向進行逐層延拓(如圖1黃色部分所示),每延拓一層,相應的延拓網格尺寸按指數增加,以保證延拓網格的稀疏性。

圖1 求解域四邊形網格剖分示意圖Fig.1 Schematic diagram of the quadrilateral meshing of the solution domain

注意,這里面的網格單元有2種比較特殊(如圖2所示)：第1種是軸單元,其所在的邊有一個落在旋轉軸上;第2種是邊界單元,其所在的邊有一個處于延拓的最外部。這2種單元的處理均與其他普通單元不同,普通單元的插值和積分均采用GLL(Gauss-Lobatto-Legendre)點。軸單元由于旋轉軸使得r→0,這在數值積分時會引入異常,因此在涉及軸單元的位置,歸一化r帶來的權重函數恰好方便采用GLJ(Gauss-Lobatto-Jacobi)積分,詳細見后文。邊界單元為了應用無限元映射法,致使將無窮遠作為插值節點,因此對應地選擇GRL積分[15],避免選擇無窮遠作為計算自由度。

圖2 單元類型以及對應的插值和積分節點選擇Fig.2 Cell type and corresponding interpolation and integration node selection

1.3 等參元映射與無限元方法

對于每一個四邊形單元,通過等參元映射的辦法可以將原四邊形單元映射到一個[-1,1]×[-1,1]的標準正方形上(如圖3(a)所示),其中位置x=(r,z)可以用等參元ξ=(ξ,η)表示,其中ξ,η∈[-1,1]。

圖3 等參元映射示意圖Fig.3 Schematic diagram of isometric mappings

這樣,對于每一個單元Ωe所示可以通過等參元定義形函數

(7)

對于普通單元(如圖3(b)所示),(ξα,ηβ)α,β=0,…,N取GLL點,這樣插值形函數Nα(ξ)可以采用Lagrange插值多項式。

而對于邊界單元,其外邊界并沒有達到無窮遠處。無限元映射法[9]可以很好地解決該問題(如圖3(c)所示),其形函數定義為

(8)

其中等參元與位置的對應關系如下

(9)

這樣,無限元映射就將無窮遠點映射到了等參元中的右端點。換言之,在等參元坐標系下,無窮遠點變成了正常的插值點,因此可以直接使用零邊界條件,然后套用譜元法求解。

1.4 插值與積分

對于任意函數f(x)可以用Lagrange插值多項式進行擬合

(10)

其偏導數可以寫為

(11)

對應插值公式(10),任意函數f(x)的積分可以寫為如下形式

(12)

其中ωi,ωj為對應的積分權重。

不同正交多項式對應的插值節點和積分權重不同。表1給出了本文用到的3種方法對應的插值節點和積分權重的值(多項式取2階)。

1.5 單元剛度矩陣計算

單元剛度矩陣可以寫成

(13)

類似地可以得到對應的單元載荷向量Fe

(14)

(15)

1.6 軸單元

由于體積元(6)中存在因子r,使得應用傳統的GLL積分在左端點為零,從而缺少了對左端點物理量的約束。為了解決該問題,Gerritsma和Phillips[14]引入了歸一化半徑的方法。例如,使用權函數ω(ξ)=1+ξ對r進行歸一化處理,方程(13)可以寫為

(16)

當ξ→-1時,η→0和ω→0,由此利用洛必達(L’Hpital’s)法則有

(17)

(18)

而剩下的ω恰好為GLJ積分的權函數。

1.7 組剛與求解

經過上述步驟后,每一個單元都可以得到一個單元剛度矩陣Ke和對應的單元載荷向量Fe。按照自由度將單元剛度矩陣Ke和單元載荷向量Fe組合起來,得到總體剛度矩陣K和總體載荷向量F。根據方程(5),K與F滿足方程

KU=F.

(19)

其中U對應要求的自由度。由于K的稀疏性,式(19)可以使用共軛梯度(conjugate gradient,CG)算法有效求解。

2 數值試驗

本文以均勻實心球(半徑為R=500 km)為例(具體參數的物理意義可以參看圖1),計算其引力位Ψ,Ψ滿足方程(1)。當rR時,ρ=0。

上述問題的解析解[16]為：

(20)

考慮到該問題滿足軸對稱性條件,可以采用上文描述的網格剖分方法(如圖1所示)。整個實心球對應的半圓(D形區,如圖1藍色部分所示)域剖分為1 630個四邊形單元。針對實心球外的真空區域,網格分別按照2、5、7、10層的層數對外部網格進行延拓(如圖1黃色部分所示),最遠延拓到約8.6倍實心球半徑處。除去最開始延拓需要對近地表網格加密,后面每向外延拓一層,網格增加60個單元,增加量不足原問題的4%。

運用上文的無限元-譜元混合法進行數值模擬,結果如圖4所示：當求解域貼近實心球表面時(藍色實線),模擬結果并不能很好地收斂于解析解(黑色虛線)。隨著外部稀疏網格層數(從2層增加到10層)的增加,數值模擬結果逐步收斂于解析解。當延拓網格達到7層(綠色實線)約3.3倍實心球半徑時,基本能滿足模擬的計算需求,進一步達到10層(紅色實線) 約8.6倍實心球半徑時,模擬結果與解析解吻合。

圖4 數值試驗模擬結果對比Fig.4 Comparison of numerical test simulation results

3 結論與討論

通過數值試驗,驗證了本文提出的無限元-譜元混合法在引力位計算模擬中的真實準確性。只需要有限的幾層網格延拓,并不需要真的把求解域延伸到無窮遠,即可實現對無窮遠邊界條件的把控。同時引入的2.5維近似假設可以很好地解決3維真實地球模型計算的復雜度。對于傳統的PREM模型[12]可以直接應用此方法。對于更復雜的地球模型(如S40RTS[11]),可以考慮保留其有效的球諧展開階次,利用傅里葉變換轉化為本文的2.5維問題[13]。最后,本文通過無限元-譜元混合法的引入,彌補了傳統引力位邊界條件近似處理的缺陷,為地球自由振蕩和低頻大尺度(如全球尺度)地震波場數值模擬提供有效支持。