基于Maxwell 模型的線性粘彈性波傳播問題的質量集中雜交應力有限元法

周亞婧, 陳豫眉

(1.四川大學數學學院， 成都 610064；2.西華師范大學公共數學學院， 南充 637009)

設Ω?R2為凸有界開集, 邊界為?Ω,T為正常數. 考慮如下基于Maxwell 模型的二維粘彈性固體介質波傳播問題：

(1)

這里u=(u1,u2)T為位移向量；ut=?u/?t,utt=?2u/?t2；σ=(σij)2×2為對稱應力張量；div為散度算子；ε(ut)=(?ut+?Tut)/2為應變張量；f=(f1,f2)T為體積力；C為彈性模量.

其中μ和λ為Lamé常數；I為二階單位矩陣；φ0，φ1和ψ0分別為位移、速度和應力的初始值. 如非特別申明，本文假設φ0，φ1，ψ0以及f的正則性是足夠的.

使用有限元法求解的過程中，時常涉及質量矩陣的求逆. 對質量矩陣進行質量集中后，得到的對角質量陣方便求逆，可以大大提高計算效率.所謂質量集中是指通過選取合適的數值求積公式將質量矩陣對角化. 眾多研究表明：對某些有限元格式可以實現質量集中[1-7]. 其中，文獻[3-4]針對矩形和長方體有限元，利用 Gauss-Lobatto 積分節點作基函數節點實現了質量集中，且其數值結果表明數值積分不影響有限元格式的數值精度.

另一方面，關于粘彈性問題的雜交應力有限元法，已有不少成果. 如文獻[8]研究了求解基于 Maxwell 模型的粘彈性波傳播問題的全離散雜交應力有限元法，分析了半離散和全離散格式的解的存在唯一性，給出相應的誤差估計；文獻[9]研究了彈性動力學問題的質量集中雜交應力有限元法，給出在一定網格條件下該數值積分在四邊形網格上的截斷誤差估計.

本文基于文獻[8]中雜交應力有限元法的半離散格式，提出了一種新的全離散方法，并研究其質量集中格式.我們在空間方向使用雜交應力四邊形有限元離散，位移采用等參雙線性插值，應力采用Pian&Sumihara的5參數應力模式，時間方向采用中心差分離散.

2 弱問題

設r為非負整數，Hr(Ω)為標準r階Sobolev空間,其范數和半范數分別記為‖·‖r和|·|r.特別地，H0=L2表示平方可積函數空間. 在不引起混淆的情況下，將與時間相關的函數v(x,t)簡記為v(t). 設k為非負整數，1≤p<∞.定義

{v:[0,T]→Hr(Ω);v關于t是k階連續可

導的}，

{v∈L1(0,T;L2(Ω));‖v‖Lp(Hr)<∞}.

當1≤p<∞時，

當p=∞時，

‖v‖L∞(Hr)=ess sup0≤t≤T‖v(·,t)‖r.

相應半范數為

及

定義如下位移空間和應力空間:

為τ=(τij)2×2的跡. Σ的范數‖·‖0定義為

下面我們給出波傳播問題式(1)的弱形式：給定初值φ0,φ1∈V,ψ0∈Σ及體積力f∈C0([0,T];(L2(Ω))2)，求(σ,u)∈C0([0,T];Σ)×C2([0,T];V)，使得

(2)

其中

(·,·)表示空間Ω上的L2內積, 即

對τ∈Σ，引入能量范數‖·‖a，即

易得

顯然，有如下的連續性條件成立：

|a(σ,τ)|‖σ‖0‖τ‖0,

?σ,τ∈C0([0,T];Σ)；

|b(τ,v)|‖τ‖0|v|1,

注1為方便起見，本文用記號ab表示存在一個和空間網格參數h，時間步長Δt均無關的正常數C使得不等式a≤Cb成立.此外，以下穩定性條件成立[10]:

(A2) Inf-sup條件：?v∈V,有

關于弱問題式(2)解的存在唯一性參見文獻[8].

3 雜交應力有限元方法

3.1 區域剖分

假設Ω是一個凸多邊形區域. 令Th表示區域Ω的一個四邊形網格剖分，其中h代表網格中所有四邊形單元直徑的最大值. 對任意四邊形單元K∈Th，其直徑記為hk. 令Zi(xi,yi),1≤i≤4表示四邊形單元K的四個頂點. 假設T1,T2,T3以及T4分別為單元K的四個子三角形Z1Z2Z3，Z1Z2Z4，Z1Z3Z4以及Z2Z3Z4.定義

假設Th滿足如下正規性條件[11]: 存在正常數ρ使得對所有的K∈Th,滿足

以“弘揚法治精神，建設法治國土”為主題，以打造“法治國土”為主線，在用好新媒體打造法治國土建設的道路上迸出法治宣傳的新火花，適應新媒體新技術發展帶來的新情況、新挑戰，著力構建全媒體普法體系，以“互聯網+（加）傳統媒體宣傳”的疊加模式，著力提升法治宣傳教育的吸引力和認可度，不斷提升慈溪市國土資源局普法依法治理水平，推動“互聯網+國土管理+全民參與+全媒體聯動”普法工作格局的建立。

hK≤ρK

(3)

(4)

對任意的四邊形單元K∈Th，定義等參雙線性映射：

x=a0+a1ξ+a2η+a12ξη,

y=b0+b1ξ+b2η+b12ξη,

其中, 單元幾何參數

(5)

注2當單元K為平行四邊形時，我們有a12=b12=0，且等參雙線性變換FK退化為一個仿射變換.

3.2 半離散/全離散雜交應力格式

令Σh，Vh分別表示應力空間Σ和位移空間V的有限維逼近，并設φ0,h，φ1,h以及ψ0,h分別為初始值φ0，φ1以及ψ0的某種逼近.弱問題式(2)的半離散格式為：

求(σh,uh)∈C0([0,T];Σh)×C2([0,T];h)滿足

(6)

對任意關于t的函數wh，記

對位移的逼近采用等參雙線性插值

?K∈Th}

(7)

K∈Th}，

其中

(8)

a1,a2,b1,b2由式(5)給出. 下述離散形式的穩定性條件成立[10]:

注3對于Stokes元Q1-P0而言，唯一不穩定情形是棋盤模式網格.后文中均假設Th為非棋盤模式網格.這意味著離散的強制性條件(B1)總是成立的.

對于半離散格式(6)，其解存在唯一，并有以下誤差估計成立[8]：

定理3.1設(σ(t),u(t))為弱問題式(2)的解，(σh(t),uh(t))為半離散格式(6)滿足初值條件

ψ0,h=IΣhψ0，φ1,h=IVhφ1

的解,這里IΣh:Σ→Σh,IVh:V→Vh分別為某種插值或投影算子，滿足

‖ψ0-IΣhψ0‖0h‖ψ0‖1,

‖φ1-IVhφ1‖0h‖φ1‖1.

若

σ∈L∞([0,T];(H1(Ω))2×2),

σt∈L2([0,T];(H1(Ω))2×2),

ut∈L∞([0,T];(H2(Ω))2),

utt∈L2([0,T];(H2(Ω))2),

則下列不等式成立：

Ch(‖φ1‖1+‖σ‖L∞(H1)+‖σt‖L2(H1)+

‖ut‖L∞(H2)+‖utt‖L2(H2)),

其中

對于以t為自變量的函數φ，定義

于是，基于二階中心差分時間離散的全離散雜交應力有限元格式為：

(9)

其中Ih表示從空間V映到Vh的等參雙線性插值算子.

證明 ?K∈Th,設

P由式(8)給出，Ni(i=1,2,…,4)為節點基函數，ξ，η是局部等參坐標，

(10)

在式(9)中的第二個等式中取任意βK∈R5，τh|K=PβK得

n≥0

(11)

引入記號

顯然HK是對稱正定矩陣.由式(11)可得

(12)

由式(9)的第一個等式可得

記

則

(13)

un=(v1,…,vm)αn,

φn=(θ1,…,θr)βn，

L(i,j)=(vi,vj),

F(i,j)=(Δt)2(vi,fn)，

Lα1=2Lα0-Lα-1-Mβ0+F.

注意，對任意的K∈Th，QK和HK都是對稱正定矩陣，從而L也是對稱正定的，而α-1,α0,β0可以由初始值ψ0,h，φ0,h得到.因此，由式(13)可以求出α1，再由式(12)可求得β1.

類似地，當2≤k≤M時，可以由(αk-1,βk-1)和(αk-2,βk-2)求出(αk,βk).定理得證.

4 質量集中雜交應力有限元法

(14)

(15)

(16)

這里JK為等參變換FK的Jacobian行列式. 由文獻[9]的結果，我們得到定理4.1.

定理4.1假設剖分 Th滿足正規性條件式(3)以及(B)條件[6]，那么對于Gauss-Lobatto積分公式(14)，下述截斷誤差估計成立:

?v,w∈Vh，?K∈Th

(17)

在全離散格式的第一個方程

?vh∈Vh

這里

同時，局部的質量矩陣

接下來，對n+1時刻的局部質量矩陣作質量集中處理. 利用Gauss-Lobatto數值積分式(15)來近似計算MK.記

將所有單元對應的局部質量集中矩陣進行組裝，所得總質量矩陣也為對角陣.

最后，全離散格式(9)的質量集中格式為：

(18)

其中

IK(·)為采用Gauss-Lobatto 公式(14)的數值積分.

注4由定理4.1和文獻[8]中使用的標準分析，我們可以證明采用了Gauss-Lobatto數值積分的質量集中格式(18)與原格式(9)的收斂階相當，且有如下誤差估計：

h+(Δt)2.

然而，數值實驗表明，由數值積分導致的相容性誤差會隨著時間步迭代逐次累積.

本節給出算例來驗證全離散雜交應力有限元法式(18)的數值性能. 取Ω=[0,1]×[0,1]，T=1.對空間區域使用圖1所示的兩種N×N四邊形網格：正方形網格和不規則四邊形網格; 對時間區間[0,T],使用M等分網格. 作為對比，我們同時給出未做質量集中的全離散雜交應力有限元格式(9)(簡記為 HSFEM)的計算結果. 為考察精度，引入下述應力和速度誤差記:

圖1 區域Ω的網格剖分Fig.1 Grid mesh Th in region Ω

例5.1設真解u=(u1,u2)T,σ=(σij)2×2具有如下形式:

u1=u2=-e-tsin(πx)sin(πy)，

σ11=3μπte-tcos(πx)sin(πy)+

λπte-tsin(πx)cos(πy)，

σ12=σ21=μπte-t(sin(πx)cos(πy)+

cos(πx)sin(πy)),

σ22=3μπte-tcos(πy)sin(πx)+

λπte-tsin(πx)cos(πy).

表1 固定時間步長Δt=0.02，μ=1，λ=1時的誤差

表2 同步加密時的誤差

表3 矩形網格μ=0.4，λ=1時的誤差