王娜，邱亞南，劉東杰

（上海大學(xué) 理學(xué)院，上海 200444）

0 引言

設(shè)Ω為具有Lipschitz 邊界Γ的單連通的有界區(qū)域，Ωc：=R2/Ω為其外部區(qū)域，對(duì)于給定的函數(shù)f∈L2(Ω)，u0∈H1/2(Γ)，t0∈H－1/2(Γ)，考察如下的界面?zhèn)鬏攩栴}：求(u，uc)∈H1(Ω) ×H1(Ωc)滿足

其中，b∈R是有界常數(shù)，n是邊界上由內(nèi)部Ω指向Ωc的法向量。本文采用Sobolev 空間中的常用記號(hào)［1］，Sobolev 空間Wm，2(Ω)簡(jiǎn)記為Hm(Ω)，||·||Hm(Ω)和|·|Hm(Ω)分別表示Hm(Ω)空間中的范數(shù)和半范數(shù)。Hm(Ω)（m是整數(shù)）的跡空間記為Hm－1/2(Γ)，跡算子γ：Hs(Ω) →Hs－1/2(Γ)，u|Γ：=γu。A?B表示存在常數(shù)C使得A≤CB。?v∈H1(Ω)，定義能量范數(shù)的形式如下：

對(duì)于上述界面問題的求解有很多不同的耦合方法［2-3］，而由馮康和余德浩教授提出的自然邊界元與有限元方法是基于相同的變分原理，可以自然而直接的耦合［4］。該耦合法的剛度矩陣就是有限元和邊界元的矩陣之和，相比有限元與古典邊界元耦合法的計(jì)算更加簡(jiǎn)單［5-6］。 采用均勻剖分的自然邊界元與有限元的耦合法已經(jīng)得到一定的應(yīng)用［7］，文獻(xiàn)［8-9］分別利用非匹配網(wǎng)格和最速下降法給出在均勻剖分下該問題的數(shù)值解。

對(duì)于有奇異性的問題，均勻剖分比較浪費(fèi)資源，而近些年比較流行的自適應(yīng)技術(shù)是處理大梯度和有間斷解問題的有效工具。自適應(yīng)技術(shù)在有限元與古典邊界元耦合中的應(yīng)用已經(jīng)比較成熟，其關(guān)鍵是找到合適的后驗(yàn)誤差估計(jì)指導(dǎo)局部網(wǎng)格細(xì)化過程，從而使數(shù)值解更快地逼近真實(shí)解，例如：

·h-h/2 后驗(yàn)誤差：此后驗(yàn)誤差估計(jì)獨(dú)立于問題，節(jié)省開銷從而備受關(guān)注。對(duì)于有限元問題，邊界元問題以及耦合問題，基于h-h/2 后驗(yàn)誤差估計(jì)的自適應(yīng)算法在飽和性假設(shè)的條件下都被證明是收斂的［10-11］。

·二水平后驗(yàn)誤差估計(jì)：Mund 等［12］利用層次基技術(shù)給出二水平后驗(yàn)誤差估計(jì)，此后驗(yàn)誤差估計(jì)方便邊界元和有限元實(shí)施獨(dú)立的細(xì)化，后來Kerb 結(jié)合Steklov-Poincare 算子將二水平后驗(yàn)誤差估計(jì)作為局部誤差指示因子來研究自適應(yīng)有限元與邊界元耦合［13］。

·基于殘差的后驗(yàn)誤差估計(jì)：Carstensen 等［14］運(yùn)用自適應(yīng)有限元與古典邊界元耦合法求解界面問題時(shí)，導(dǎo)出了一個(gè)基于殘差的后驗(yàn)誤差估計(jì)給出了誤差的一個(gè)上界。

對(duì)于線性界面?zhèn)鬏攩栴}，h-h/2 后驗(yàn)誤差估計(jì)等價(jià)于二水平后驗(yàn)誤差估計(jì)［11］。Feischl 等［15］針對(duì)界面?zhèn)鬏攩栴}分別討論基于二水平后驗(yàn)誤差估計(jì)和基于殘差的后驗(yàn)誤差估計(jì)并分析了自適應(yīng)過程的收斂率，文獻(xiàn)［16］給出兩種后驗(yàn)誤差估計(jì)的可靠性和有效性分析。

綜上，關(guān)于自適應(yīng)有限元與古典邊界元耦合已有較多的研究成果，而關(guān)于自適應(yīng)自然邊界元的耦合研究較少，本文利用自然邊界元的良好性質(zhì)，研究自適應(yīng)有限元與自然邊界元耦合法。為了敘述和計(jì)算的方便，我們考慮Ω為圓域的線性界面?zhèn)鬏攩栴}（1）的自適應(yīng)有限元與自然邊界元耦合法。具體框架如下：第一部分給出了等價(jià)的變分問題；第二部分給出離散的變分問題；第三部分給出本文的兩個(gè)誤差分析，即h-h/2后驗(yàn)誤差估計(jì)和基于殘差的后驗(yàn)誤差估計(jì)；最后利用數(shù)值算例驗(yàn)證我們的理論結(jié)果。

1 變分問題

針對(duì)半徑為R的圓外區(qū)域問題，引入自然積分算子K：H1/2(Γ) →H－1/2(Γ)和泊松積分算子P：H1/2(Γ) →H1(Ω)［5］，則泊松積分方程和自然積分方程如下所示，* 號(hào)是指自然積分算子K 作用于函數(shù)u0，

設(shè)(u，uc)是方程（1）的解，對(duì)方程（1）的第一個(gè)表達(dá)式兩邊同乘v∈H1(Ω)，應(yīng)用Green 公式可得

將方程（1）的第5 個(gè)表達(dá)式代入（5）式得

由方程（1）的第二個(gè)式子和第四個(gè)式子可得如下的自然積分方程和泊松積分方程，

利用自然積分算子K 的雙線性，將（7）式代入（6）中得到，

也即(u，uc) ∈H1(Ω) ×H1(Ωc)滿足如下的變分問題，

因此，變分問題（10）等價(jià)于微分方程（1）。顯然，求解變分問題時(shí)，只需求出u∈H1(Ω)即可，而uc依賴于u|Γ利用（10）的第二個(gè)表達(dá)式可以求出。

我們定義算子A(u，v) ：H1(Ω) ×H1(Ω) →R和線性泛函L(v) ：H1(Ω) →R的表達(dá)式如下：

則變分問題（10）等價(jià)于下述方程，

由商空間中的等價(jià)模定理可知［1］，存在常數(shù)C，使得對(duì)?u∈H1(Ω)/P0，

因此，雙線性型A(u，u) 也是V-橢圓的，存在常數(shù)C使得

又因?yàn)殡p線性型A(u，u) 是連續(xù)的［4］，由Lax-Milgram 定理可得，變分問題（13）在商空間H1(Ω) /P0中存在唯一解。

2 離散變分問題

設(shè)Th是區(qū)域Ω的擬正則三角剖分，Th/2是對(duì)Th的一致細(xì)化得到的網(wǎng)格剖分，Th+1是對(duì)Th實(shí)施一次自適應(yīng)加密細(xì)化后得到的新一層網(wǎng)格剖分。Eh表示剖分三角單元所有邊的集合，Eh(Ω)和Eh(Γ)分別為內(nèi)部和邊界上的邊的集合。N 記為剖分節(jié)點(diǎn)集合，N (Ω)和N (Γ)分別表示區(qū)域內(nèi)和邊界上所有節(jié)點(diǎn)的集合，n1，n2分別為集合N (Ω)，N (Γ)中元素的個(gè)數(shù)。對(duì)于任意的三角單元T∈Th，hT∶= diam(T)，Eh(T)是單元T的三邊的集合，對(duì)任意的邊界單元E∈Eh，hE∶= diam(E)，h?！? maxE∈Eh(hE)，h?！謍E≈hT≈|T|1/2。

為了便于實(shí)施，令邊界的剖分與區(qū)域的三角剖分在邊界上是一致的。原問題（13）的離散變分形式如下：

3 后驗(yàn)誤差估計(jì)

后驗(yàn)誤差主要利用數(shù)值解和已知條件估計(jì)真實(shí)解和數(shù)值解的誤差大小。對(duì)于我們研究的線性問題，h－h(huán)/2 后驗(yàn)誤差估計(jì)和二水平后驗(yàn)誤差估計(jì)是等價(jià)的［10］，因此本節(jié)主要研究h－h(huán)/2 法和殘差法這兩個(gè)主流的后驗(yàn)誤差估計(jì)。

3.1 h－h(huán)/2后驗(yàn)誤差估計(jì)

本節(jié)主要在有限元與自然邊界元的耦合框架下利用h－h(huán)/2 后驗(yàn)誤差估計(jì)的思想給出與誤差的能量范數(shù)|||u－uh|||等價(jià)的估計(jì)量ηl和μl

令數(shù)值解滿足飽和性假設(shè)條件：

由（17）式和（18）式易得，存在常數(shù)C使得后驗(yàn)誤差估計(jì)量ηl和誤差|||u－uh||| 滿足下述關(guān)系

存在常數(shù)C∞＞0 使得?T∈Th，

那么利用上面兩個(gè)式子可得

其中，T′ ∈Th/2是一致細(xì)化Th中的三角單元T得到的任一子單元，且有|T| ≤v|T′|，v為常數(shù)。

定理1 在飽和性假設(shè)的條件下，存在常數(shù)C1，C2＞0，使得

由Glerkin 解的擬最優(yōu)性可知

由文獻(xiàn)［18］引理2.2 可得，

由最后兩個(gè)估計(jì)式（24）式和（25）式得ηl?μl。

綜上可得，μl?ηl?μl，結(jié)合（19）式可得

因此，μl是誤差|||u－uh|||的有效且可靠的后驗(yàn)誤差估計(jì)量。

3.2 基于殘差的后驗(yàn)誤差估計(jì)

本節(jié)主要利用自然積分算子的性質(zhì)，推導(dǎo)出一個(gè)形式較為簡(jiǎn)單且易于計(jì)算的后驗(yàn)誤差估計(jì)。

其中ηh(T)的表達(dá)式如下：

令[?uh·n]表示uh在單元的邊界E上法向?qū)?shù)的躍度，對(duì)區(qū)域Ω每個(gè)元素T∈Th上應(yīng)用Green 公式可得，

將（32）式代入（31）式可得

由Ho¨lder 不 等 式 和 引 理1 可 得，

再利用Cauchy 不等式可得

結(jié)合引理1 和引理2 可得對(duì)任意三角單元T的邊E∈Eh有下述不等式成立，

從而可得（33）式右端第二項(xiàng)和第三項(xiàng)有如下估計(jì)，

將以上不等式（35）－（37）代入（33）式可得

4 數(shù)值實(shí)驗(yàn)

本節(jié)主要利用上面給出的兩個(gè)后驗(yàn)誤差估計(jì)進(jìn)行自適應(yīng)有限元與自然邊界元耦合的數(shù)值實(shí)驗(yàn)。自適應(yīng)流程如下：

輸入：初始網(wǎng)格，自適應(yīng)細(xì)化參數(shù)θ∈(0，1)

（1）求出相應(yīng)的變分問題的數(shù)值解uh；

（2）計(jì)算的每個(gè)網(wǎng)格單元上的后驗(yàn)誤差估計(jì)(μl(T)2+μl(E)2)或ηh(T)；

（3）將誤差大的單元τ標(biāo)記組成集合M，使得

（4）細(xì)化被標(biāo)記的元素得到新的網(wǎng)格剖分；

重復(fù)這個(gè)過程：求解→估計(jì)→標(biāo)記→細(xì)化，直到數(shù)值解的誤差或者細(xì)化網(wǎng)格的直徑足夠小。

4.1 數(shù)值解法

4.2 誤差范數(shù)的計(jì)算

設(shè)u和uh是方程（13）和（14）的解，我們記LHS（Th1）和LHS（Th2）分別為定理1 和定理2 的誤差項(xiàng)，GUB（Th1）和GUB（Th2）分別記為定理1 和定理2 中的后驗(yàn)誤差估計(jì)，

在計(jì)算誤差LHS（Th1）和LHS（Th2）時(shí)，由于范數(shù)‖·‖H1/2(Ω)難以計(jì)算，將||u|Γ－uh|Γ||H1/2(Γ)利用其等價(jià)的范數(shù)‖hΓ(u|Γ－uh|Γ)′‖L2(Γ)近似，而‖u－uh‖H1(Ω)和GUB（Th1）和GUB（Th2）可利用數(shù)值求積公式計(jì)算。

4.3 實(shí)驗(yàn)結(jié)果

我們?nèi)ˇ笧閱挝粓A，精確解u和uc的取為（41）式，其滿足界面?zhèn)鬏敺匠蹋?）。初始網(wǎng)格剖分取為Ω的內(nèi)接正六邊形，將其平均等分為六個(gè)正三角元素，在邊界Γ上用直線近似，自適應(yīng)細(xì)化參數(shù)取為θ=0.5。

圖1（a）給出了基于h－h(huán)/2 后驗(yàn)誤差估計(jì)的自適應(yīng)加密網(wǎng)格的三角剖分圖，從圖上可以看出自適應(yīng)網(wǎng)格加密主要集中在真實(shí)解u波動(dòng)較大的第三象限的區(qū)域。圖1（b）給出了隨著網(wǎng)格一致加密或自適應(yīng)加密真實(shí)誤差和后驗(yàn)誤差估計(jì)量的收斂結(jié)果，從圖上可以看出一致加密和自適應(yīng)的結(jié)果都是收斂的，且自適應(yīng)的收斂速度比一致細(xì)化得更快，說明我們的算法是有效的。

圖1 基于h－h(huán)/2 后驗(yàn)誤差估計(jì)的自適應(yīng)結(jié)果（a）基于GUB（Th1）的自適應(yīng)網(wǎng)格剖分圖；（b）基于GUB（Th1）的自適應(yīng)誤差收斂圖Fig. 1 Adaptive results of the h－h(huán)/2 posteriori error estimation（a）Adaptive grid generation diagram based on GUB（Th1）；（b）Error convergence graph based on GUB（Th1）

從圖2（a）自適應(yīng)剖分的結(jié)果可以看出，自適應(yīng)網(wǎng)格加密主要集中在真實(shí)解波動(dòng)較大的第三象限的區(qū)域和部分邊界，這也符合我們的預(yù)期結(jié)果。圖2（b）是一致剖分和自適應(yīng)算法的誤差收斂結(jié)果，從圖上可以看出隨著網(wǎng)格的細(xì)化，數(shù)值解和真實(shí)解的誤差線性的收斂到零，且自適應(yīng)誤差下降速率比一致細(xì)化的快。

圖2 基于殘差的后驗(yàn)誤差估計(jì)的自適應(yīng)結(jié)果（a）基于GUB（Th2）的自適應(yīng)網(wǎng)格剖分圖；（b）基于GUB（Th2）的自適應(yīng)誤差收斂圖Fig. 2 Adaptive results of a posteriori error estimation based on residuals（a）Adaptive grid generation diagram based on GUB（Th2）；（b）Error convergence graph based on GUB（Th2）

從上述兩個(gè)后驗(yàn)誤差的數(shù)值實(shí)驗(yàn)結(jié)果來看，h－h(huán)/2 后驗(yàn)誤差估計(jì)的計(jì)算比較簡(jiǎn)單，而基于殘差的后驗(yàn)誤差估計(jì)使誤差收斂地更快。