虛擬元方法求解混合邊界條件下的線彈性問題

馬俊馳, 程新博, 梁曉坤

(遼寧師范大學 數學學院,遼寧 大連 116081)

線彈性問題在材料科學、橋梁結構、土木工程和航空航天[1-3]等領域具有廣泛的應用.目前線彈性問題已有多種數值求解方法,如有限元方法[4]、有限差分法[5]、有限體積法[6]、弱有限元方法[7]、虛擬元方法[8]等.有限元方法在我國最早是由馮康先生[9]提出的,并形成一套系統的數值計算方法.有限元方法概念清晰、描述簡單、應用廣泛,但在剖分單元的選擇上具有一定的局限性,比如多邊形(多面體)網格,而虛擬元方法[10]在網格處理上具有較好的靈活性,其作為有限元方法的推廣受到廣泛關注.本文主要研究應用虛擬元方法求解二維混合形式線彈性問題,對于該問題可以加強應力張量的對稱性,這種方法的挑戰在于混合形式線彈性問題的離散化.文獻[11]中對應力張量的對稱性進行了弱化,引入旋轉張量作為另一個未知數,主要討論具有齊次Dirichlet邊界條件的二維線彈性問題.而本文將邊界條件推廣到Dirichlet和Neumann的混合邊界條件.文獻[12]研究具有混合邊界條件的二維線彈性問題的三種弱形式,給出虛擬元離散化的一般過程,并沒有給出詳細的虛擬元空間構造、誤差估計等,本文將研究具有混合邊界條件的弱對稱形式的線彈性問題.

1 連續性問題

定義Ω?2為具有Lipschitz連續邊界的有界區域,邊界為Γ=?Ω,且Γ的兩個子集ΓD,ΓN滿足Γ=ΓD∪ΓN,ΓD∩ΓN=?.記彈性體所受外力的單位體密度為f∈[L2(Ω)]2.線彈性問題可表述為,求σ,u滿足

(1)

2 虛擬元方法

2.1 變分形式

定義空間Σ:=H(div;Ω),U:=[L2(Ω)]2,X:={Υ:Υ∈[L2(Ω)]2×2,Υ反對稱}.

通過分部積分法,得到問題(1)的變分形式:找到(σ,u,ω)∈Σ×U×X,使得

(2)

據文獻[13]得雙線性形式a(·,·)滿足連續性和強制性,b(·,·),c(·,·)滿足inf-sup條件,則問題(2)的解存在且唯一.

2.2 構造虛擬元空間

首先對區域Ω進行剖分,設{Th}h是將Ω分解為元素E的序列,h表示{Th}h中元素的最大直徑.假設存在整數N和正實數β,使得對于每個h和每個E∈{Th}h,滿足:①E的邊數小于等于N;②E的最短邊與E的直徑hE之比大于β;③E相對于半徑為βhE的球上的每個點呈星形.

對于整數k≥1,定義局部應力空間為

為了對誤差進行估計,引入下列投影算子與插值算子.

(3)

(4)

對于足夠光滑的函數q和l,以下估計成立

(5)

在每個單元E上,

(6)

2.3 虛擬元離散

基于上述剖分以及空間的定義,首先將雙線性形式進行離散,

由于雙線性形式b(·,·),c(·,·)可以通過選擇適當的空間Σh,Uh,Xh來計算,所以不需要引入全局項b(·,·),c(·,·)的任何近似值,即bh(·,·)=b(·,·),ch(·,·)=c(·,·).

變分形式(2)的離散格式為:找到(σh,uh,ωh)∈Σh×Uh×Xh,使得

(7)

其中,div Σh?Uh.

上述的局部離散雙線性形式具有以下的3種性質.

引理1[16]存在整數k≥1,對于所有的E∈Th,都有[Pk(E)]2×2?Σh|E,

(2) 穩定性:存在兩個不依賴于h和E的正常數α*和α*,使得

根據引理可知離散雙線性形式ah(·,·)具有連續性、強制性,b(·,·),c(·,·)滿足離散inf-sup條件,則離散形式(7)的解存在且唯一.

2.4 誤差分析

根據文獻[11]可得離散問題的收斂性,下面對離散問題的誤差進行分析.

引理2[13]存在一個常數C>0,使得

對于誤差分析,給出如下定理并證明.

定理1問題(7)有唯一解(σh,uh,ωh)∈Σh×Uh×Xh,L2誤差估計為

‖σ-σh‖0+‖u-uh‖0+‖ω-ωh‖0≤C(μ)(hk+1‖σ‖k+1+hk+1‖u‖k+1+hk+1‖ω‖k+1),

其中,C(μ)是一個與λ無關,但依賴于μ的正常數.

根據問題(2)和問題(7),為了方便起見,這里取gD=0,有

ah(eσ,τh)+b(τh,eu)+c(τh,eω)=

(8)

(9)

(10)

(11)

再根據式(11)有

(12)

令引理2中vh=0,θh=eω得到τh∈Σh滿足c(τh,eω)≥C‖τh‖0‖eω‖0,由式(8)有

(13)

根據式(12),式(13)得到估計

接下來,根據三角不等式和式(5),式(6)可以推得

‖σ-σh‖0+‖ω-ωh‖0≤C(μ)(hk+1‖σ‖k+1+hk+1‖ω‖k+1).

(14)

(15)

根據問題(15)以及式(4),式(8)和τ的對稱性有

其中,

接下來由一致性、連續性和式(14)可得

同理可得

R2≤C(μ)(hk+2‖σ‖k+1+hk+2‖ω‖k+1)‖eu‖0,

R3≤C(μ)(hk+2‖σ‖k+1+hk+2‖ω‖k+1)‖eu‖0,

將R1,R2,R3,R4代回原式,兩邊同時除以‖eu‖0就可以得到

(16)

由三角不等式和式(16)可得

‖u-uh‖0≤C(μ)(hk+1‖σ‖k+1+hk+1‖u‖k+1+hk+1‖ω‖k+1).

結論得證.

3 數值實驗

本文考慮帶孔的板模型,其中,孔的半徑為a=0.4.考慮混合弱形式線彈性問題(1),其中,彈性模量為E=103,泊松比為υ=0.3.弱對稱形式經過剛度矩陣的變換可化為強對稱形式,即以位移和應力作為問題的解,因此數值實驗中只考慮位移與應力.根據文獻[17]可得,位移的精確解為

應力的精確解為

其中,(r,θ)為極坐標,θ是從x軸的正半軸開始沿逆時針方向的夾角.對于平面應變狀態,κ=3-4υ.μ為Lamé常數.

首先對區域進行剖分,剖分區域為Ω=Ω1/Ω2,其中,

圖1給出剖分區域Ω的多邊形網格剖分,其中,網格剖分數n分別為46、539、3592和21916.

圖1 不同網格數對應的剖分圖Fig.1 The division diagram corresponding to different mesh numbers

圖2給出網格剖分數n為46、539、3592和21916的位移與應力的數值解,其中,u1,h,u2,h為位移的數值解,σ11,h為應力的數值解.

圖2 不同網格數對應的位移與應力的數值解Fig.2 Numerical solution of displacement and stress corresponding to different mesh numbers

其次應用虛擬元方法求解上述線彈性問題,得到數值解uh和σh與精確解u和σ之間的L2范數下的相對誤差,討論不同網格剖分數下位移與應力的相對誤差.

表1給出了不同網格剖分數對應的位移與應力的相對誤差,其中,uerr(L2)表示L2范數下位移的相對誤差,σerr(L2)表示L2范數下應力的相對誤差,其中,uerr(L2)=‖u-uh‖0,σerr(L2)=‖σ-σh‖0.

表1 不同網格數對應的位移與應力的相對誤差Table 1 Relative error of displacement and stress corresponding to different mesh number

從表1可以看出,隨著網格剖分數的增加,位移與應力的相對誤差都呈明顯減小的趨勢,這與第三節的理論分析相一致,驗證了誤差估計的有效性.

4 總 結

本文主要討論帶Dirichlet和Neumann混合邊界條件的二維線彈性問題的虛擬元方法.線彈性問題研究的是在一定假設條件下,彈性體受到外力或者其他外界作用時,彈性體內的應力、應變和位移之間的關系.本文研究混合弱形式的線彈性問題,以位移、應力和旋轉張量作為研究問題的解,在混合邊界條件約束下,通過變分原理得到弱形式線彈性問題的變分格式,然后通過構造虛擬元空間和自由度對變分形式進行離散.進一步通過定義投影算子與插值算子處理誤差估計,最后通過實際問題帶孔的板,驗證了理論分析的有效性.