雙曲型界面問(wèn)題的變網(wǎng)格低階有限元方法

關(guān)宏波,洪亞鵬,曲雙紅

( 鄭州輕工業(yè)大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院,河南 鄭州450002)

1.引言

界面問(wèn)題在材料科學(xué)、流體力學(xué)和電磁場(chǎng)問(wèn)題等工程領(lǐng)域有著重要應(yīng)用,例如帶有不同傳導(dǎo)系數(shù)的熱傳導(dǎo)問(wèn)題、具有不同粘性系數(shù)的兩相流問(wèn)題以及電磁場(chǎng)中存在絕緣分界體等具有兩種或兩種以上介質(zhì)的情形[1].此類問(wèn)題的精確解通常不易解析表達(dá),于是研究其數(shù)值方法就變得十分有意義.而有限元方法在眾多數(shù)值計(jì)算方法中具有舉足輕重的地位,界面問(wèn)題的有限元方法根據(jù)網(wǎng)格剖分形式大體可以分為兩類：擬合有限元方法和非擬合有限元方法.其中非擬合有限元方法是指網(wǎng)格剖分與界面所處位置無(wú)關(guān),如具有代表性的浸入有限元方法[2],其主要技巧是修改界面單元上的基函數(shù)使其滿足界面條件,得到界面單元內(nèi)部的分片多項(xiàng)式空間,新的求解空間對(duì)界面問(wèn)題的解有較好的逼近效果[3?5].但是該方法也有不足之處,比如即使原問(wèn)題滿足相應(yīng)的連續(xù)性和強(qiáng)制橢圓性,得到的數(shù)值解卻也是非對(duì)稱的,界面單元所對(duì)應(yīng)的剛度矩陣與非界面單元?jiǎng)偠染仃嚥町愝^大.而擬合有限元方法是指在網(wǎng)格剖分時(shí)要考慮界面所處位置進(jìn)行適當(dāng)剖分,文[6]最早引入了對(duì)求解區(qū)域沿界面進(jìn)行特殊剖分,使界面作為單元邊界或單元頂點(diǎn),此類方法的優(yōu)勢(shì)在于最終合成的總剛矩陣與傳統(tǒng)偏微分方程結(jié)果一樣,均為對(duì)稱正定的,這給數(shù)值計(jì)算帶來(lái)諸多便利之處.基于此,文[7]提出了一種無(wú)限元方法,用無(wú)限剖分的思想與有限元方法的結(jié)合,得到了能量模意義下的最優(yōu)誤差估計(jì)結(jié)果,但是該方法對(duì)于曲邊型界面問(wèn)題是不適用的.文[8]克服以上缺陷,得到了有限元方法的最優(yōu)誤差估計(jì)結(jié)果,但是其對(duì)精確解的正則性假設(shè)過(guò)高,不但要求精確解及其外法向?qū)?shù)沿界面方向連續(xù),而且要求解在每個(gè)子區(qū)域里具有四階偏導(dǎo)數(shù),這一限制過(guò)于苛刻.陳志明院士和鄒軍教授在文[9]中提出了一個(gè)簡(jiǎn)單的P1協(xié)調(diào)三角形有限元方法,在L2模和能量模意義下得到了次優(yōu)的誤差估計(jì)結(jié)果,并且在文中的一個(gè)注解里給出了能夠達(dá)到最優(yōu)誤差估計(jì)的一個(gè)正則性假設(shè)條件.鄒軍教授等人[10]在對(duì)此類問(wèn)題解的正則性進(jìn)行更為精細(xì)探討的基礎(chǔ)上,通過(guò)引入一個(gè)擬插值算子,得到了一般k階有限元逼近下的最優(yōu)階誤差估計(jì)結(jié)果.本文作者之一將其推廣到了非協(xié)調(diào)有限元逼近格式[11],同樣導(dǎo)出了最優(yōu)階誤差估計(jì).兩個(gè)印度學(xué)者Sinha和Deka對(duì)文[9]中關(guān)于假設(shè)正則性得到更高收斂階的注解進(jìn)行了重新論證,得到了最優(yōu)階L2模和能量模誤差估計(jì),并應(yīng)用于拋物型界面問(wèn)題的半離散和全離散Galerkin格式[12?13].然而在許多實(shí)際問(wèn)題中,界面的位置和形狀都有可能會(huì)隨著時(shí)間的不同而產(chǎn)生相應(yīng)的變化,于是在不同時(shí)間節(jié)點(diǎn)處允許采用合適的網(wǎng)格剖分具有一定的應(yīng)用價(jià)值.梁國(guó)平院士在文[14]中針對(duì)拋物問(wèn)題提出了一種變網(wǎng)格有限元處理方法,其基本做法是:對(duì)空間域采用有限元方法,對(duì)時(shí)間軸采用差分法,并且對(duì)于不同時(shí)間的空間區(qū)域可以采用不同的網(wǎng)格.后續(xù)有不少國(guó)內(nèi)外學(xué)者將該方法應(yīng)用到更多的非定常偏微分方程問(wèn)題中,其中包括協(xié)調(diào)和非協(xié)調(diào)有限元方法[15?18].然而對(duì)于界面問(wèn)題的變網(wǎng)格方法,目前尚未見(jiàn)到有文獻(xiàn)正式出版.本文將采用變網(wǎng)格技術(shù),將擬合有限元方法應(yīng)用到雙曲型界面問(wèn)題中,采用最低階的P1三角形協(xié)調(diào)有限元對(duì)空間變量進(jìn)行逼近,得到相對(duì)于網(wǎng)格剖分尺寸和時(shí)間離散步長(zhǎng)的最優(yōu)誤差估計(jì)結(jié)果,并將該方法進(jìn)行推廣應(yīng)用于其他情形,從而拓寬了有限元的應(yīng)用范圍.

2.雙曲型界面問(wèn)題

考慮如下雙曲型界面問(wèn)題:求u ∈L2([0,T];X(?)),使得

其中? ?R2為一個(gè)凸多邊形區(qū)域,??為?內(nèi)部的一個(gè)開(kāi)區(qū)域,其邊界Γ=???也包含于?內(nèi),并且具有C2光滑性,?+=????,如下圖2.1所示.定義X(?)=H10(?)∩H2(?+)∩H2(??),本文中所用到的Sobolev空間及模均為標(biāo)準(zhǔn)記法(見(jiàn)文[19]).

圖2.1 ?=??∪?+

問(wèn)題(2.1)在界面處滿足如下跳躍性條件:

其中[u]Γ為u跨過(guò)界面Γ的跳度,n為相應(yīng)于界面Γ的單元外法向量,而

這里s=+或?.

令Q=,(2.1)式等價(jià)于下面形式:

(2.4)式所對(duì)應(yīng)的變分形式為：求u(x,t)∈(0,T]→H10(?),使得

其中

下面對(duì)區(qū)域?進(jìn)行擬一致三角剖分Th.首先用一個(gè)凸多邊形??h逼近??,??h的頂點(diǎn)都位于Γ上,?+h=????h.然后對(duì)??h和?+h分別進(jìn)行擬一致三角剖分,但是其中的三角形單元K要滿足以下兩個(gè)條件：

1)K要么在??h內(nèi),要么在?+h內(nèi);

2)假設(shè)F為K的一條邊,如果F ∩?！?,則F的兩個(gè)頂點(diǎn)在Γ上或者整條邊都在Γ上.

如果K與Γ沒(méi)有交點(diǎn),稱此類單元為非界面單元; 否則K為一個(gè)界面單元,記界面單元的全體為Th?,而hK為單元K的直徑,對(duì)于任意的界面單元K,記K?=K ∩??,K+=K ∩?+.由于Γ是C2光滑的,則meas(K?)≤ch3K與meas(K+)≤ch3K中必然有一個(gè)是成立的,其中meas(·)為測(cè)度,這里及以后出現(xiàn)的c均為一個(gè)與網(wǎng)格剖分尺寸及時(shí)間離散步長(zhǎng)無(wú)關(guān)的正常數(shù).為方便起見(jiàn),我們記K為Ks中滿足meas(Ks)≤ch3K的那一個(gè).

在Th上定義線性協(xié)調(diào)P1三角形有限元空間Vh,并記Πh:H10(?)→Vh為相應(yīng)的插值算子,由文[9]知,如果進(jìn)一步假設(shè)?0為包含界面的某小一鄰域,當(dāng)h足夠小時(shí),必然有T?h ??0,滿足u ∈W1,∞(??∩?0)∩W1,∞(?+∩?0),能夠得到如下最優(yōu)階的插值誤差估計(jì):

(2.5)式的有限元離散形式為:求uh,Qh ∈Vh滿足

3.變網(wǎng)格方法及誤差估計(jì)

首先引入變網(wǎng)格的主要思想,設(shè)0=t0< t1< t2< ··· < tN=T將時(shí)間軸平均分成N段,分點(diǎn)為tn,n=0,1,2,··· ,N,?t=tn+1?tn=令Thn為tn時(shí)刻對(duì)空間區(qū)域?的一個(gè)三角形剖分簇,V nh為tn時(shí)刻的有限元空間:V nh= {v(x,tn);v ∈Vh}.

我們選取uh(x,t)的近似解空間Sh如下:Sh上的函數(shù)uh(x,t),有N+1個(gè)有限元插值函數(shù),即uh(x,t)在N+1個(gè)節(jié)點(diǎn)的節(jié)點(diǎn)值,在時(shí)間間隔tn < t < tn+1內(nèi),uh(x,t)取為t的線性函數(shù),通過(guò)線性連結(jié)相鄰兩個(gè)節(jié)點(diǎn)uh(x,tn)和uh(x,tn+1)而得到.近似解uh(x,t)在各個(gè)時(shí)間點(diǎn)處的函數(shù)值unh=uh(x,tn)由以下格式確定:

變網(wǎng)格有限元解與精確解的誤差主要包含以下三個(gè)部分:對(duì)時(shí)間變量的差分誤差、對(duì)空間變量的有限元誤差以及網(wǎng)格變動(dòng)誤差.下面的引理3.1首先給出差分誤差的估計(jì).

引理3.1記un=u(x,tn),Qn=Q(x,tn),則成立

證注意到Vh ?H10(?)為協(xié)調(diào)有限元空間,在(2.5)中令v=vh,并從tn到tn+1積分,得

然后(3.2)與(3.4)相減,可得

利用差分性質(zhì)及Cauchy不等式,上式右端兩項(xiàng)可分別估計(jì)為:

和

引理得證.

為書寫方便,我們引入以下記號(hào):

其中Πnun為tn時(shí)刻un的有限元插值,n=1,2,··· ,N?1.

引理3.2假設(shè)時(shí)間剖分參數(shù)0

其中

證(3.1)中第五個(gè)式子與(3.2)相減,有

由(3.7)得

從而結(jié)(3.10),有

在上式中取vh=ξn+1+,則

由(3.1)中最后一個(gè)式子得到

將上式代入(3.12)式得

其中β?=min(β+,β?).

我們對(duì)上式右端五項(xiàng)利用Cauchy不等式,擬不等式及引理3.1,分別估計(jì)為:

和

將(3.14)-(3.18)結(jié)果代入(3.13)式中,有

由(3.1)式中第一個(gè)式子,顯然(?(?rn),?vh)=(?(?en),?vh),?vh ∈V hn+1.在上式中取vh=,則||1≤|rn|1+|?en|1,對(duì)上式兩端平方,并由Cauchy不等式可得

同樣利用上述技巧,由(3.1)式中第二個(gè)式子可得到

將(3.20)和(3.21)式代入(3.19)式,有

而當(dāng)兩層網(wǎng)格之間不發(fā)生變動(dòng)時(shí),為

假設(shè)計(jì)算過(guò)程中共有M次網(wǎng)格變動(dòng),不妨設(shè)前M層網(wǎng)格變動(dòng),后N?M層網(wǎng)格不變動(dòng),將各層對(duì)應(yīng)的估計(jì)式寫出,并適當(dāng)使用技巧對(duì)其求和,可得

下面給出本文主要結(jié)論:

定理3.1假設(shè)unh,Qnh ∈V nh和un,Qn ∈X(?)分別為(2.5)和(2.7)的解,則

證因?yàn)?/p>

所以由引理3.2得

引理得證.

注3.1由估計(jì)式可以看出,網(wǎng)格的變動(dòng)次數(shù)不能是任意的,只有當(dāng)M為一個(gè)與h及?t無(wú)關(guān)的有界數(shù)時(shí),誤差的估計(jì)階才是最優(yōu)的.

4.推廣應(yīng)用

本文主要針對(duì)雙曲型界面問(wèn)題提出了一個(gè)變網(wǎng)格的有限元方法,在理論分析及推廣應(yīng)用過(guò)程中,需要指出以下三點(diǎn):

1)在誤差估計(jì)過(guò)程中,我們用到了該P(yáng)1有限元空間的關(guān)鍵性質(zhì),即對(duì)于任意的vh ∈Vh,(?vh)|K為常數(shù).于是該方法對(duì)于非協(xié)調(diào)的P1三角形元[11]和CNQrot1矩形元[20]也是適用的,但對(duì)于其它一些常見(jiàn)的協(xié)調(diào)和非協(xié)調(diào)元,該理論還不能直接應(yīng)用.

2)該方法可以通過(guò)適當(dāng)調(diào)整逼近格式,推廣應(yīng)用于其它的一些界面問(wèn)題,如半線性界面問(wèn)題、線彈性界面問(wèn)題等.

3)當(dāng)定義區(qū)域?為一個(gè)窄邊型區(qū)域時(shí),可以拋棄擬一致剖分限制,考慮采用各向異性剖分[21],從而達(dá)到減少總體自由度的效果.

總之,該變網(wǎng)格方法為處理界面問(wèn)題提供了一種新的思路,對(duì)設(shè)計(jì)時(shí)間相關(guān)界面問(wèn)題的自適應(yīng)算法是有益的.