朱雙彪
(南京財(cái)經(jīng)大學(xué)應(yīng)用數(shù)學(xué)學(xué)院,江蘇 南京 210023)
許多科學(xué)問題都可以轉(zhuǎn)化為偏微分方程問題.通過數(shù)值求解這類問題的方法有很多,自然邊界元[1-3]與有限元耦合法就是其中的一類.耦合法的好處在于,2種方法結(jié)合時能夠彌補(bǔ)各自的不足.利用耦合法不僅可以求解圓和橢圓區(qū)域問題,還可以求解凹角區(qū)域問題[4-5].但在求解耦合問題時,由于是用直邊三角形作有限元分析的,因此會產(chǎn)生誤差.采用曲邊有限元即曲邊三角形代替直邊三角形[6-8],就可以大大降低誤差.筆者擬給出凹角區(qū)域上求解泊松方程邊值問題的曲邊有限元和自然邊界元的耦合法.
如圖1所示,取夾角α、線段Γ0和Γ1、光滑的圓弧Γ,一起圍成區(qū)域Ω.在區(qū)域Ω內(nèi)求解如下邊值問題:
圖1 區(qū)域ΩFig. 1 Area Ω
(1)
(2)
用半徑為R的人工圓弧Γ2置于區(qū)域Ω內(nèi),其中Γ2={(R,θ)|0≤θ≤α},且dist(Γ,Γ2)>0.Γ2將區(qū)域Ω分為2個區(qū)域,Γ2包圍的區(qū)域?yàn)棣?,Γ2和Γ包圍的區(qū)域?yàn)棣?,如圖2所示.在區(qū)域Ω2內(nèi)用自然邊界元方法,在區(qū)域Ω1內(nèi)用曲邊有限元方法.
圖2 區(qū)域Ω引入人工邊界Fig. 2 Artificial Boundary of Area Ω
在耦合法中,一般是用直邊三角形近似劃分有限區(qū)域,但這種處理方式存在一定的誤差.采用曲邊元素,就可以很好地?cái)M合待求解區(qū)域彎曲的邊界[4-5].
(3)
取參考單元G,其頂點(diǎn)g1=(0,0),g2=(1,0),g3=(0,1).記面積坐標(biāo)(λ1,λ2,λ3),于是得到如下線性變換:
其中λ1=1-ξ-η,λ2=ξ,λ3=η.作如下變換F:
(4)
其中
(5)
(6)
這里
線性變換(3)~(6)可以將參考單元的3條邊一一映射為曲邊三角形的2條直邊和1個曲邊[6].
σ={Ni:i=1,2,3},是一組線性函數(shù).定義Ni(φ)=φ(pi),φ∈Σ,i=1,2,3.在τ和有限元正則的情況下,各個曲邊三角單元的插值誤差估計(jì)為
(7)
變分問題(2)的離散問題是
(8)
變分問題(8)有唯一的連續(xù)依賴于f的解.
引理1[1]
證明因?yàn)镾h∈H1(Ω1),所以在
中取v=vh∈Sh,再與
相減,得到
又因?yàn)閡h∈Sh,所以
證畢.
由引理1和(7)式可以得到定理1,2:
定理1(收斂性) 假設(shè)插值算子滿足
那么
定理2(能量模估計(jì)) 設(shè)u∈H2(Ω1),那么
由定理1給出的近似解滿足收斂性可知,通過自然邊界元和曲邊有限元耦合法求解問題(1)是可行性.
‖u-uh‖L2(Ω1)≤Ch2‖u‖2,Ω1.
對比能量模估計(jì)和L2模估計(jì),L2模估計(jì)是最優(yōu)的.
例1在圖1中取夾角α=7π/4.采用曲邊有限元與自然邊界元耦合法求解如下外邊值問題:
其中Γ0,Γ1,Γ的邊界條件分別為
Γ0={(r,θ)|r≥3,θ=0},Γ1={(r,θ)|r≥3,θ=α},Γ={(r,θ)|r=3,0<θ<α},
人工邊界為
Γ2={(r,θ)|r=1.5,0<θ<α}.
通過Matlab軟件對區(qū)域進(jìn)行加密剖分,結(jié)果如圖3,4所示.圖3是將邊界劃分為33個節(jié)點(diǎn),圖4是將邊界劃分為17個節(jié)點(diǎn),可以看出圖3的剖分比圖4更密.
圖3 mesh=8×33Fig. 3 mesh=8×33
圖4 mesh=4×17Fig. 4 mesh=4×17
取剖分網(wǎng)格N=17,33,65,129.有限元采用直邊三角元和曲邊三角元產(chǎn)生的誤差分別見表1,2.
表1 有限元采用直邊三角元產(chǎn)生的誤差Table 1 Error Value of Straight Triangular Element in Finite Element
表2 有限元采用曲邊三角元產(chǎn)生的誤差Table 2 Error Value of Curved Edge Element in Finite Element
從表1,2可知:誤差隨著剖分網(wǎng)格的加密越來越??;當(dāng)剖分網(wǎng)格相同時,L2誤差階都接近于ο(h2),但有限元采用曲邊三角元產(chǎn)生的誤差比直邊三角元的小.這說明,曲邊有限元采用曲邊三角單元代替原來的直邊三角形,可以克服邊界近似產(chǎn)生的誤差,降低數(shù)值積分誤差.