付燦燦, 紀海峰

(南京郵電大學(xué)理學(xué)院, 南京 210000)

界面問題是在研究不同材料或相同材料不同狀態(tài)之間相互作用過程中提出的, 如研究復(fù)合材料中的熱傳導(dǎo)問題[1-3]、彈性界面問題[4]、多相流問題[2,5]等．根據(jù)物理守恒定律,界面問題需要在交界面上滿足一些跳躍條件, 而不同問題所要滿足的界面條件各不相同．眾所周知, 除了某些簡單的方程外, 偏微分方程的解析求解一般都比較困難[6],因此如何設(shè)計快速、有效的數(shù)值求解方法對解決實際問題有重要意義．然而,對于界面問題,其復(fù)雜的界面跳躍條件給數(shù)值方法的設(shè)計帶來了不小的挑戰(zhàn)．目前求解界面問題的數(shù)值方法按網(wǎng)格剖分大致可以歸為兩類:匹配網(wǎng)格方法和非匹配網(wǎng)格方法．所謂匹配網(wǎng)格,是指網(wǎng)格必須界面互相匹配, 即界面只能穿過單元的頂點,不能穿過單元內(nèi)部．對于這種匹配的網(wǎng)格, 有限元方法能夠捕捉真解在界面上的信息,可獲得最優(yōu)逼近性[7-9]．然而,對復(fù)雜界面, 尤其是三維情形, 生成界面匹配網(wǎng)格非常困難[2], 且當(dāng)界面隨時間移動時在求解過程中該方法需要在每一時間層上重新生成界面匹配網(wǎng)格,這不僅增加了計算量,而且給理論分析帶來困難．為了克服上述缺點,有學(xué)者提出非匹配網(wǎng)格方法,即界面與網(wǎng)格無關(guān),且允許其以任意方式穿過單元內(nèi)部．為了獲得最優(yōu)逼近性,非匹配網(wǎng)格方法必須在界面單元(即界面穿過的單元)上進行精細處理[10-13]．浸入有限元方法是一類簡單有效的非匹配網(wǎng)格方法[14-17],它通過修改界面單元上的基函數(shù)獲得最優(yōu)的收斂性,但這種修改帶來了一個負面影響,即基函數(shù)在界面邊(界面穿過的邊)上不連續(xù),這給有限元方法帶來了相容誤差．為了消除該誤差, Zhuang等[16-17]提出了部分罰浸入有限元方法．該方法利用間斷有限元思想,在界面邊上添加額外的懲罰項來懲罰基函數(shù)的不連續(xù)性,而理論分析指出,罰參數(shù)需要充分大才能保證該方法的穩(wěn)定性,且在實際計算時罰參數(shù)需要人為猜測選取,這給實際計算帶來麻煩．為克服罰參數(shù)選取的缺陷, 針對二階橢圓界面問題, Ji等[18]提出了參數(shù)友好型部分罰浸入有限元方法．該方法通過構(gòu)造提升算子并把其添加到懲罰項中來保證穩(wěn)定性,其優(yōu)點是不需要人為選擇罰參數(shù)即可保證算法的穩(wěn)定性．受其啟發(fā),本文擬把上述參數(shù)友好型部分罰浸入有限元方法推廣到拋物界面問題并進行理論分析．在空間方向上將采用參數(shù)友好型部分罰浸入有限元方法離散, 而在時間方向上則采用經(jīng)典的Crank-Nicolson格式, 并借助橢圓投影算子, 推導(dǎo)半離散及全離散格式的最優(yōu)誤差估計．

1 拋物界面問題

設(shè)Ω是R2上的凸多邊形域,Γ是浸入Ω中的C2光滑界面, 且把Ω分成2個不相交的子區(qū)域Ω+和Ω-．不失一般性, 假設(shè)Ω-嚴格包含在Ω的內(nèi)部, 如圖1所示．考慮如下拋物界面問題:

(1)

圖1 區(qū)域Ω=Ω+Ω-ΓFig.1 Domain Ω=Ω+Ω-Γ

(2)

u(·,t)=0, 在?Ω上,t>0,

(3)

其解滿足跳躍條件

(4)

其中n為Γ的單位法向量并指向Ω+; [·]Γ表示在界面Γ上的跳躍, 即[v]Γ=(v+-v-)|Γ, 式中v+=v|Ω+,v-=v|Ω-; 系數(shù)β在Ω上為分片常數(shù),即

2 浸入有限元方法

設(shè){Τh}h>0是Ω上的一族三角剖分, 對于每一個單元T∈Th, 定義其網(wǎng)格參數(shù)為h=max{hT|T∈Th}, 其中hT為單元T的直徑．假設(shè)網(wǎng)格滿足以下條件:

(H1) 界面Γ與單元的邊交點數(shù)目不超過2, 除非這條邊是界面Γ的一部分．

(H2) 若界面Γ和單元有2個交點, 那么這2個交點一定在這個單元2條不同的邊上．

(H3) 剖分Τh是正則的．

圖2 界面單元Fig.2 Interface element

在界面單元上的形函數(shù)定義為:

(5)

且滿足

(6)

其中nh是Γh的單位法向量．令Sh(T)是形函數(shù)在界面單元T上組成的空間, 由定義可知, 形函數(shù)φ∈Sh(T)及其通量在Γh∩T上連續(xù), 即

(7)

若界面單元的最大角αmax≤π/2, 則函數(shù)φ∈Sh(T)由φ(Ai)(i=1,2,3)唯一確定[18], 其中Ai(i=1,2,3)為T的頂點．

(8)

則式(1)～(3)的半離散方程為

(9)

(10)

其中

(11)

(12)

(13)

(14)

(15)

(16)

3 誤差分析

引理1[18]存在正常數(shù)C1,C2, 使

Ah(u,v)≤C1|‖u|‖h·|‖v|‖h,

(17)

(18)

Ah(Rhv,ωh)=Ah(v,ωh),?ωh∈Sh．

(19)

(20)

(21)

(22)

4 例證

為驗證收斂階, 構(gòu)造數(shù)值算例: 考慮區(qū)域Ω=(-1,1)×(-1,1), 界面Γ為半徑是0.5, 圓心在(0,0)處的圓, 真解為

圖3 誤差收斂示意圖Fig.3 Error convergence diagram

其中β+=10,β-=1．

首先將區(qū)域分成N×N個矩形, 再把每個矩形沿對角線分成兩個三角形, 從而得到一致三角剖分且網(wǎng)格尺寸為2/N．取T=1, 并在時間方向進行等距剖分, 步長選為Δt=h．圖3給出了誤差收斂示意圖．從圖3可以看出,L2誤差與斜率為2的直線基本一致．根據(jù)定理1, 在T時刻有數(shù)值解的L2誤差應(yīng)為2階,驗證了本文所提方法的正確性．而H1誤差與斜率為1的直線基本一致,故H1誤差為1階收斂．