董 成,鄔吉明

(1.中國工程物理研究院研究生院,北京 100088)

(2.北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所,北京 100088)

1 引言

考慮擴(kuò)散方程邊值問題

其中??R2是平面上的多邊形區(qū)域,Λ是擴(kuò)散張量,f是源項(xiàng),??=ˉΓD∪ˉΓN是?的邊界,n為單位外法向量,gD和gN分別是Dirichlet和Neumann邊值.

上述擴(kuò)散問題有著廣泛的實(shí)際應(yīng)用背景,例如在輻射流體力學(xué)、輻射磁流體力學(xué)、油藏模擬等應(yīng)用中,擴(kuò)散過程與其他物理過程相耦合.我們常常需要在網(wǎng)格強(qiáng)扭曲、強(qiáng)間斷、強(qiáng)各項(xiàng)異性、強(qiáng)非線性等極端條件下求解擴(kuò)散方程,這是一件具有挑戰(zhàn)性的工作.有限體積方法是求解擴(kuò)散問題最常用的方法之一,它具有局部守恒、簡單容易實(shí)現(xiàn)等良好性質(zhì).多年來,許多科研工作者致力于擴(kuò)散方程的有限體積方法研究,提出了眾多的離散格式,如九點(diǎn)格式(NPS)[1]、多點(diǎn)通量逼近格式(MPFA)[2]、支撐算子格式[3]、廣義差分格式[4]等.按照未知量的類型,這些格式可大致分為單元中心格式、節(jié)點(diǎn)格式、雜交格式、混合格式等等,詳細(xì)的最新研究進(jìn)展可參見文獻(xiàn)[5–9].

我國學(xué)者李德元教授等人在二維網(wǎng)格上基于積分插值法構(gòu)造了一個(gè)單元中心型有限體積格式[1],由于該格式在結(jié)構(gòu)四邊形網(wǎng)格上有九點(diǎn)模板,故常稱其為九點(diǎn)格式,它是若干輻射流體力學(xué)程序的基本格式[10,11].由于單元中心格式在一個(gè)單元上只有一個(gè)未知量,在構(gòu)造離散格式時(shí)需要引入輔助未知量來提高精度.九點(diǎn)格式的輔助未知量定義在網(wǎng)格節(jié)點(diǎn)處,如何用單元中心未知量逼近節(jié)點(diǎn)輔助未知量是九點(diǎn)格式研究中的一個(gè)重要內(nèi)容.一個(gè)理想的二維節(jié)點(diǎn)插值算法應(yīng)當(dāng)具有簡單、保正、不依賴于網(wǎng)格拓?fù)?、不依賴于間斷線的位置、有二階精度、易于向三維推廣等性質(zhì).目前已知的節(jié)點(diǎn)插值算法,如算術(shù)平均插值[1]、逆距離加權(quán)插值[12]、最小二乘插值[13]、顯式加權(quán)算法[14]等,均只能滿足部分要求.由于缺乏容易向三維推廣的二階插值算法,嚴(yán)重制約了九點(diǎn)格式在三維問題中的應(yīng)用.

本文在多點(diǎn)通量逼近的邊輔助未知量插值算法的基礎(chǔ)上,通過一個(gè)特殊的極限技巧獲得了一個(gè)新的節(jié)點(diǎn)插值算法,并在給定假設(shè)下證明了該算法中局部線性系統(tǒng)的可解性.本文插值算法的一個(gè)重要之處在于容易向三維推廣,且滿足除保正性以外的其他要求.全文的推導(dǎo)基于所謂的線性精確準(zhǔn)則[15],即當(dāng)擴(kuò)散張量關(guān)于網(wǎng)格是分片常數(shù)、解析解關(guān)于網(wǎng)格是分片線性函數(shù)時(shí),算法推導(dǎo)的每一步都是精確成立的.為了和通常的等式相區(qū)別,用符號’表示相關(guān)等式是在線性精確意義下成立的.

2 九點(diǎn)格式的構(gòu)造

最初的九點(diǎn)格式是在擴(kuò)散系數(shù)為標(biāo)量的情形下推導(dǎo)的,但我們?nèi)菀讓⑵渫茝V到擴(kuò)散系數(shù)為張量的情形[14],為保持本文內(nèi)容的完整性,我們簡要介紹一下擴(kuò)散張量情形下九點(diǎn)格式的一種構(gòu)造算法．首先將?剖分為互不重疊的多邊形網(wǎng)格,K表示其中的一個(gè)單元,σ表示單元邊,K 的所有邊組成的集合記為EK,單元K 關(guān)于邊σ的單位外法向量記為nK,σ.記流向量為F=?Λ?u,這里假定Λ關(guān)于網(wǎng)格是分片常數(shù),并用ΛK表示Λ在單元K 上的限制.單元K 通過邊σ的流記為FK,σ,即

在單元K上對方程(1.1)兩端積分并利用散度定理,有

于是有限體積法解擴(kuò)散問題就歸結(jié)為構(gòu)造流FK,σ的離散表達(dá)式.

在圖1中,K和L表示兩個(gè)相鄰單元,它們的單元中心分別是OK和OL,σ是它們的公共邊,其頂點(diǎn)為A 和B,nK,σ,nL,σ分別是K,L關(guān)于σ的單位外法向量.dK,σ,dL,σ分別表示OK,OL到σ的距離,uA,uB,uK,uL表示u在點(diǎn)A,B,OK,OL處的近似值.對于向量v=(a,b)T,記v?=(b,?a)T.此外,為敘述簡單起見,引入以下記號(見圖1)

首先,有如下的向量分解

其中

圖1:離散流構(gòu)造模板

|σ|表示σ的長度.其次,在單元K 上,根據(jù)文獻(xiàn)[14]的引理2.1,

注意到 |σ|nK,σ=t?,t=aK? bK,·t= ?|σ|dK,σ,將 (2.4) 和 (2.5) 式代入 (2.1) 式,經(jīng)過簡單的計(jì)算后就可以得到

類似地,在單元L上,有

通常要求通過邊σ的流滿足連續(xù)條件,即

由此得到

將(2.6)和(2.7)式代入(2.9)式的右端,利用aK?bK=aL?bL=t,aK?aL=bK?bL=s,并經(jīng)過簡單化簡就得到內(nèi)部邊σ上流的離散表達(dá)式

其中

容易看出,上述離散流涉及單元中心未知量uK和uL以及節(jié)點(diǎn)輔助未知量uA和uB,要得到純單元中心格式,需要用單元中心量對節(jié)點(diǎn)輔助未知量進(jìn)行插值.

3 節(jié)點(diǎn)輔助未知量的插值算法

如圖2(a),假設(shè)Q0是一內(nèi)部節(jié)點(diǎn),它周圍的單元?k(1≤k≤n)按逆時(shí)針順序排列.Q0Pk和Q0Pk+1是?k的以Q0為頂點(diǎn)的邊,Ok為其中心.假設(shè)Tk是邊Q0Pk上的一點(diǎn),由下式確定

其中τ∈(0,1).令u0,uk,分別表示Q0,Ok,Tk處的未知量,Λk表示Λ在?k上所取的分片常數(shù)值,nk表示?k關(guān)于Q0Pk的單位外法向,見圖2(b).本文中k將被視為以n為周期的指標(biāo),例如當(dāng)k=n+1和0時(shí),有Pn+1=P1,T0=Tn.

圖2:內(nèi)部節(jié)點(diǎn)Q0處的局部結(jié)構(gòu)和記號

3.1 邊輔助未知量的插值算法

在多點(diǎn)通量逼近格式中,邊未知量ˉuk可以通過單元中心未知量uk表示.為敘述簡潔起見,引入以下記號

其中ek表示n階單位矩陣的第k個(gè)列向量.顯然,Tk是一個(gè)n×2的矩陣,并且

在單元?k的邊Q0Pk和Q0Pk+1上,定義如下離散流

注意到對于 Q0Pk(Q0Pk+1),有于是

其中

這里記號(τ)表示相關(guān)的量是τ的函數(shù),在不引起歧義的情況下將其省略.利用(3.3)式,可將(3.4)式改寫為

再利用邊Q0Pk+1上的流連續(xù)條件,即

或者其等價(jià)形式

就得到

其中

通過求解(3.8)式,就能夠用單元中心未知量來表示邊未知量,即=M?1NU.容易看出,只有當(dāng)M可逆時(shí),上述邊未知量插值算法才是可行的.自多點(diǎn)通量逼近算法提出來的二十多年里,人們還沒有在數(shù)值計(jì)算中遇到過M奇異的情形,但這一點(diǎn)至今仍缺乏嚴(yán)格的理論證明.

3.2 新的節(jié)點(diǎn)未知量插值算法

其中Sk,1(Sk,2)表示4Q0PkOk(4Q0OkPk+1)的面積.于是(3.5)式可以寫為

其中

再利用(3.9)式,可得

其中

需要注意,由于(3.11)式的分母中含有τ,Ml和Nl依然是τ的函數(shù).根據(jù)(3.11)式,發(fā)現(xiàn)

故有

其中1和0分別表示所有分量為1和0的向量,O表示零矩陣.利用(3.12)和(3.14)式,可證(3.8)式等價(jià)于

注意到解沿著Q0Pk(1≤k≤n)是切向連續(xù)可微的,通過泰勒展開,可以得到

其中 Γ =(γ1,···,γn)T是與 τ無關(guān)的常向量,O(τ2)表示所有分量為 O(τ2)的向量. 將(3.16)式代入(3.15)式并且利用(3.14)式中的第一個(gè)等式,得到

最后將以上方程兩端同時(shí)除以τ,并令τ→0,就得到最終的方程

其中M1(0),M2(0)和N2(0)表示M1,M2和N2取τ=0的值.將M1(0)的第一列替換為M2(0)1其余列不變得到一個(gè)新的矩陣,記為,則(3.17)式可以寫為與τ無關(guān)的線性系統(tǒng)

求解上述局部線性系統(tǒng)就得到一個(gè)新的節(jié)點(diǎn)插值算法,即

這里需要特別強(qiáng)調(diào),雖然上述算法基于多點(diǎn)通量逼近的邊輔助未知量插值,但它與(3.8)式中矩陣M的可逆性沒有關(guān)系.

4 局部線性系統(tǒng)(3.18)的可解性

雖然(3.8)式的可解性目前尚無理論證明,但局部線性系統(tǒng)(3.18)式的可解性在一定條件下是可以嚴(yán)格證明的.本節(jié)的主要結(jié)果依賴以下兩個(gè)假設(shè):

在擴(kuò)散系數(shù)是標(biāo)量的情況下,假設(shè)(H1)和(H2)的幾何意義是明顯的,此時(shí),兩個(gè)假設(shè)變?yōu)?/p>

圖3:假設(shè)(H1)和(H2)的幾何意義

引理4.1在假設(shè)(H1)下,M2(0)1的所有元素都是負(fù)的,即

證 由M2(0)的定義,有

其中δij表示Kronecker delta函數(shù).利用(3.11)式并且注意到τ=0,進(jìn)一步得到

再結(jié)合假設(shè)(H1),立即得到(4.1)式.

引理4.2記A=(aij)為n×n的矩陣,滿足

其中ak,bk(1≤k≤n)是正數(shù).則

證 根據(jù)(4.2)式,有

分裂An?1的最后一列,可得

對于An?1的第一部分,通過將最后一列加到第(n?2)列,再將新的第(n?2)列加到第(n?3)列,依此類推,最終可得到一個(gè)上三角型行列式,進(jìn)而有

用同樣的方法可以得到

注意到A1=a1+bn>0以及ai,bi>0,通過數(shù)學(xué)歸納法最終得到An?1>0,從而>0.類似地,可以證明>0.當(dāng)i 6=1,n時(shí),記aii對應(yīng)的余子式為,并令

其中Ik是k階單位矩陣.容易驗(yàn)證是形如(4.4)式的三對角矩陣.于是

定理4.3 在假設(shè)(H1)和(H2)的條件下,(3.18)式中的系數(shù)矩陣 eM是非奇異的.

證 由M1(0)的定義和(3.13)式,有

由(3.11)式和假設(shè)(H2),可以得到

令A(yù)=(aij)n×n= ?M1(0).根據(jù)引理4.2,>0(1≤k≤n).由 eM 的定義,并利用引理4.1,可得

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

首先,定義如下的L2誤差和流誤差

其中M為單元集合,eK表示單元K中心的誤差,Nσ表示邊σ的相鄰單元數(shù),L表示σ的除K 之外的另一相鄰單元（如果存在的話）,當(dāng)σ在邊界上時(shí),我們約定dL,σ=eL=0.兩個(gè)網(wǎng)格層之間的收斂率Rα(α=u,q)通過以下公式獲得

其中h1,h2代表兩個(gè)網(wǎng)格層的網(wǎng)格尺寸,Eα(h1),Eα(h2)為對應(yīng)的離散誤差.

用BICGSTAB方法求解離散線性系統(tǒng),并取收斂準(zhǔn)則為εlin=1.0E-15.離散流采用公式(2.10),節(jié)點(diǎn)插值算法分別采用文獻(xiàn)[14]中的顯式加權(quán)算法2和本文的極限加權(quán)算法,對應(yīng)的格式簡記為NPS-EW2和NPS-LW.采用的計(jì)算網(wǎng)格見圖4,且所有的計(jì)算均在雙精度下完成.

圖4:計(jì)算網(wǎng)格

5.1 實(shí)驗(yàn)1:各向異性問題

考慮?=[0,1]2上的具有全Dirichlet邊界的擴(kuò)散問題(1.1).擴(kuò)散張量和精確解如下

這個(gè)數(shù)值實(shí)驗(yàn)是文獻(xiàn)[5]中的一個(gè)測試的一個(gè)簡單變形.

從實(shí)驗(yàn)結(jié)果來看,NPS-LW在各網(wǎng)格上與NPS-EW2的實(shí)驗(yàn)結(jié)果接近,這里僅給出前者的計(jì)算結(jié)果.表1給出了NPS-LW在各個(gè)網(wǎng)格上數(shù)值解和流的誤差,圖5以對數(shù)折線圖的形式展示了收斂速度.可以看出NPS-LW在這些網(wǎng)格上的數(shù)值解的誤差隨著網(wǎng)格加密以趨近2階的速度收斂.

表1:NPS-LW在實(shí)驗(yàn)1中的數(shù)值結(jié)果

圖5:NPS-LW在實(shí)驗(yàn)1中的誤差曲線

5.2 實(shí)驗(yàn)2:泊松問題

在[0,1]2上求解泊松方程??u=4.0,采用Dirichlet邊界條件,解析解取為

本算例十分簡單,其設(shè)計(jì)目的是測試新的節(jié)點(diǎn)加權(quán)算法與已有顯式加權(quán)算法在數(shù)值表現(xiàn)上的差異.NPS-LW的數(shù)值結(jié)果在Mesh1,Mesh2和Mesh4上與NPS-EW2非常接近,但在Mesh3上明顯優(yōu)于NPS-EW2,在表2中進(jìn)行了對比.

表2:NPS-LW與NPS-EW2在Mesh3上的數(shù)值結(jié)果對比

5.3 實(shí)驗(yàn)3:間斷系數(shù)問題

在?=[0,1]2上求解擴(kuò)散問題(1.1)–(1.2),擴(kuò)散張量取為

I為二階單位矩陣,精確解取為

這里采用Mesh1和Mesh4.對于Mesh1,所有位于x=0.5上的節(jié)點(diǎn)在x方向上都不擾動(dòng).數(shù)值計(jì)算結(jié)果見表3.從該表中可以看到,數(shù)值解誤差收斂速度均趨于2階,流的誤差接近1階,均為最優(yōu)階,表明新的節(jié)點(diǎn)加權(quán)算法也能很好地適應(yīng)間斷系數(shù)問題.

表3:實(shí)驗(yàn)3 NPS-LW在Mesh1,4上的數(shù)值結(jié)果