一類傳輸問題的自適應(yīng)FEM-BEM方法

王娜，邱亞南，劉東杰

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

0 引言

設(shè)Ω為具有Lipschitz 邊界Γ的單連通的有界區(qū)域，Ωc：=R2/Ω為其外部區(qū)域，對于給定的函數(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 空間中的常用記號［1］，Sobolev 空間Wm，2(Ω)簡記為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ù)的形式如下：

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

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

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

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

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

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

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

1 變分問題

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

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

將方程（1）的第5 個表達式代入（5）式得

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

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

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

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

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

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

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

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

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

2 離散變分問題

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

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

3 后驗誤差估計

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

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

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

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

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

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

那么利用上面兩個式子可得

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

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

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

由文獻［18］引理2.2 可得，

由最后兩個估計式（24）式和（25）式得ηl?μl。

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

因此，μl是誤差|||u－uh|||的有效且可靠的后驗誤差估計量。

3.2 基于殘差的后驗誤差估計

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

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

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

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

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

再利用Cauchy 不等式可得

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

從而可得（33）式右端第二項和第三項有如下估計，

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

4 數(shù)值實驗

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

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

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

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

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

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

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

4.1 數(shù)值解法

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

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

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

4.3 實驗結(jié)果

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

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

圖1 基于h－h(huán)/2 后驗誤差估計的自適應(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)格加密主要集中在真實解波動較大的第三象限的區(qū)域和部分邊界，這也符合我們的預(yù)期結(jié)果。圖2（b）是一致剖分和自適應(yīng)算法的誤差收斂結(jié)果，從圖上可以看出隨著網(wǎng)格的細化，數(shù)值解和真實解的誤差線性的收斂到零，且自適應(yīng)誤差下降速率比一致細化的快。

圖2 基于殘差的后驗誤差估計的自適應(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）

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