董榮榮 張 超 張耀明

(山東理工大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,山東淄博 255049)

引言

科學(xué)與工程領(lǐng)域中許多問(wèn)題，例如，穩(wěn)定熱傳導(dǎo)、彈性桿件的扭轉(zhuǎn)、穩(wěn)定滲流、動(dòng)水壓力、薄膜平衡、Helmhotz 方程、電磁場(chǎng)、非均質(zhì)材料及非線性問(wèn)題等[1-3]，的數(shù)值分析可直接地或者間接地歸結(jié)為L(zhǎng)aplace 或Poisson 方程的邊值問(wèn)題的求解.邊界元法是求解此類問(wèn)題的強(qiáng)有力的數(shù)值工具.

物理量梯度邊界積分方程由物理量邊界積分方程關(guān)于空間變量取導(dǎo)數(shù)獲得[4]，它對(duì)于某些問(wèn)題，如裂紋問(wèn)題、波散射問(wèn)題、薄板彎曲問(wèn)題及具有退化邊界的狹窄和薄域問(wèn)題等[5-7]，是非常有用的，因?yàn)榇藭r(shí)僅由基本物理量的邊界積分方程無(wú)法準(zhǔn)確地表示原邊值問(wèn)題的解，即與原邊值問(wèn)題不等價(jià).為了避免出現(xiàn)這種情況，通常將基本物理量的梯度邊界積分方程與物理量邊界積分方程組合，即對(duì)偶邊界積分方程，來(lái)表示原邊值問(wèn)題的解.位勢(shì)問(wèn)題的位勢(shì)梯度邊界積分方程已收到了許多研究[8-12].文獻(xiàn)[8]建立了二維位勢(shì)問(wèn)題的直接變量規(guī)則化梯度邊界積分方程，文獻(xiàn)[9]建立了三維位勢(shì)問(wèn)題的直接變量規(guī)則化梯度邊界積分方程.兩個(gè)梯度方程中，超奇異積分的規(guī)則化形式是相似的，即

需指出的是，u,k(y)既不是已知量，也非未知量，數(shù)值實(shí)施時(shí)需將u,k(y)沿曲面(曲線)的兩個(gè)(一個(gè))切線方向及法線方向進(jìn)行分解，切線方向的量采用插值逼近，法線方向的量作為未知量，因此這個(gè)積分在高價(jià)曲面單元下的數(shù)值實(shí)施是非常復(fù)雜和困難的.文獻(xiàn)[10]研究了二維位勢(shì)問(wèn)題的位勢(shì)梯度邊界積分方程.引入一系列變換將梯度邊界積分方程轉(zhuǎn)換為具有強(qiáng)奇異積分的自然邊界積分方程，然后采用加減技術(shù)規(guī)則化強(qiáng)奇異積分.方法有效地消去了奇異性，取得了很好的數(shù)值結(jié)果，但很難推廣到三維位勢(shì)問(wèn)題.文獻(xiàn)[11-12]建立了三維位勢(shì)問(wèn)題的間接變量規(guī)則化梯度邊界積分方程，給出了高階單元下的數(shù)值實(shí)施方案，取得了準(zhǔn)確的數(shù)值結(jié)果.方程中規(guī)則化積分形式是

它比直接法中的形式簡(jiǎn)單得多，數(shù)值實(shí)施也容易得多.但是，其數(shù)值實(shí)施也需要許多技術(shù)性措施[11-14].不同于前述的規(guī)則化邊界積分方程法，文獻(xiàn)[15-17]和文獻(xiàn)[18-21]分別提出了直接計(jì)算梯度方程中的奇異積分的局部規(guī)則化方法.它們的優(yōu)點(diǎn)是算法具有一般性，可計(jì)算任何階的奇異積分.缺點(diǎn)是兩種算法的理論推導(dǎo)是繁復(fù)的，計(jì)算量也相當(dāng)大，因?yàn)樗鼈冃鑼⒎e分核中的每個(gè)依賴于單元參數(shù)的函數(shù)表示成在內(nèi)蘊(yùn)坐標(biāo)系統(tǒng)下的局部距離ρ 的冪級(jí)數(shù).此外，還有一些別的直接計(jì)算方法[22-26]，都屬于局部法，涉及與局部單元參數(shù)有關(guān)的操作.

本文提出了求解三維位勢(shì)問(wèn)題的位勢(shì)梯度邊界積分方程的新算法.該方法通過(guò)構(gòu)造輔助邊值問(wèn)題來(lái)計(jì)算梯度邊界積分方程中的系統(tǒng)矩陣，算法實(shí)施中不需要建立規(guī)則化邊界積分方程也不需要直接計(jì)算強(qiáng)奇異積分.因此，方法具有數(shù)學(xué)理論簡(jiǎn)單、計(jì)算效率高、結(jié)果精度高等優(yōu)點(diǎn).需要強(qiáng)調(diào)的是，本文輔助邊值問(wèn)題法為梯度邊界積分方程中的強(qiáng)奇異積分計(jì)算提供了一個(gè)新的思路和途徑.另外，本文方法可以拓展到其他問(wèn)題，如彈性力學(xué)問(wèn)題、Stokes 及Helmholtz 問(wèn)題等.

1 問(wèn)題描述

本文設(shè)Ω 是示R3中的一個(gè)有界區(qū)域，Ωc是它的補(bǔ)域，Γ 是它們的邊界.n=(n1,n2)是區(qū)域Ω 的邊界Γ 在x點(diǎn)處的單位外法向量.三維位勢(shì)問(wèn)題的基本解是

具有邊界條件

這里Γ=Γ1∪Γ2且Γ1∩Γ2=φ;是邊界Γ 上已知的函數(shù).

引理[27-30]設(shè)ψ(x)∈C0,α(Γ)和是曲面Γ 上的任一光滑點(diǎn)，當(dāng)時(shí)，令，和假設(shè)(K1是一個(gè)常數(shù))，那么有

2 位勢(shì)梯度邊界積分方程的輔助邊值問(wèn)題法

2.1 邊界積分方程

現(xiàn)在從式(6)出發(fā)，導(dǎo)出位勢(shì)梯度邊界積分方程.對(duì)任一個(gè)給定點(diǎn)，由式(6)可得

2.2 八節(jié)點(diǎn)四邊形二次單元

邊界元法的實(shí)施包括邊界幾何的描述和邊界函數(shù)的插值逼近.為了一般性，本文采用八節(jié)點(diǎn)四邊形二次單元來(lái)描述邊界幾何，八節(jié)點(diǎn)二次不連續(xù)插值函數(shù)來(lái)逼近單元上的邊界函數(shù).

八節(jié)點(diǎn)四邊形二次單元的形函數(shù)Nk(ξ1,ξ2)(k=1,2,···,8)是

這里ξ1,ξ2是無(wú)因次坐標(biāo)，且-1 ≤ξ1,ξ2≤1.因此單元上的點(diǎn)x可以近似地表示為

這里xk是第k個(gè)插值節(jié)點(diǎn)的坐標(biāo).

單元上的邊界量由八節(jié)點(diǎn)二次不連續(xù)插值函數(shù)來(lái)逼近，即

這里φk是邊界函數(shù)在第k個(gè)節(jié)點(diǎn)處的函數(shù)值，α ∈(0,1)的常數(shù)，本文取參數(shù)α=2/3.

2.3 精確單元

許多情況下，邊界幾何具有參數(shù)表示.此時(shí)采用精確單元計(jì)算，可減少誤差.精確單元的概念最早是本文作者2004 年提出的[27]，后來(lái)許多學(xué)者跟隨了這個(gè)方法，甚至賦予一個(gè)新的名字，但本質(zhì)是一樣的.

在三維空間中，多數(shù)情況下，固體模型的表面可以表示成參數(shù)形式

這里a,b,c,d都是有限數(shù).當(dāng)區(qū)間[a,b]和[c,d]分別被離散成N1和M1個(gè)子區(qū)間后，相應(yīng)的曲面被離散成N×M個(gè)所謂的單元.由于這樣的分割是在參數(shù)空間內(nèi)，對(duì)應(yīng)的幾何點(diǎn)仍然在原來(lái)的曲面上，因此稱為精確單元.在每個(gè)精確單元上，參數(shù)θ,φ 可表示為

這里θi,φi(i=1,2)是精確單元的端點(diǎn)坐標(biāo).

2.4 邊界積分方程的離散化

將邊界Γ 離散成Ne單元，因而有8 ×Ne個(gè)邊界節(jié)點(diǎn)：xi，i=1,2,···,8 ×Ne.方程(7)和(9)中的y取為任一場(chǎng)點(diǎn)xi∈Γα(屬于第α 個(gè)單元)，那么方程(7)和(9)的離散化形式是

2.5 輔助邊值問(wèn)題法

數(shù)值實(shí)施中，確定式(16)中的Ai是一個(gè)困難的問(wèn)題.Ai對(duì)應(yīng)離散系統(tǒng)矩陣中的對(duì)角元，并且對(duì)角占優(yōu)，因此Ai的準(zhǔn)確與否對(duì)解系統(tǒng)的性能及解的精度影響特別大.另一方面，估計(jì)Ai需要計(jì)算一個(gè)強(qiáng)奇異積分，其強(qiáng)奇異核函數(shù)是基本解關(guān)于坐標(biāo)變量的偏導(dǎo)數(shù)，一般不是邊界法向?qū)?shù)，因此它的計(jì)算是相當(dāng)困難的[11-12].采用規(guī)則化邊界積分方程可避免直接計(jì)算此類積分[11-14]，但規(guī)則方程的形式較復(fù)雜，程序設(shè)計(jì)難度較大；直接計(jì)算此類積分需要繁復(fù)的理論推導(dǎo)和計(jì)算工作量[15-21].本文提出計(jì)算式(15),式(16)對(duì)應(yīng)的系統(tǒng)矩陣[Gjl]和[Qjl]的輔助邊值問(wèn)題法.構(gòu)造輔助超定邊值問(wèn)題，求解此邊值問(wèn)題可求得[Gjl]和[Qjl].下面給出具體實(shí)施過(guò)程:

這里(a,b,c)是Ωc外的任一點(diǎn).

現(xiàn)在用式(15)和式(16)求解邊值問(wèn)題(17)或(19)，就可獲得[Gjl]和[Qjl].值得注意的是，在求出這兩個(gè)矩陣后，使用式(15)和式(16)求解有限域Ω 或者無(wú)限域Ωc上的任何邊值問(wèn)題，不再需要計(jì)算[Gjl]和[Qjl].由此可看出，輔助邊值問(wèn)題法的效率并不差.

3 數(shù)值算例

考慮3 個(gè)數(shù)值算例，來(lái)驗(yàn)證方法的有效性.數(shù)值實(shí)驗(yàn)的重點(diǎn)在于考察輔助邊值問(wèn)題法計(jì)算邊界通量?u/?x1的能力及準(zhǔn)確性.為了估計(jì)數(shù)值誤差，采用如下L2范數(shù)

其中,M表示計(jì)算點(diǎn)數(shù)，分別是第k個(gè)計(jì)算點(diǎn)處的精確解和數(shù)值解.REu，REq及REq1分別表示邊界位勢(shì)u、法向梯度q=?u/?n及熱流通量q1=?u/?x1的平均相對(duì)誤差.

算例1作為第一個(gè)算例，考慮立方體區(qū)域Ω={(x1,x2,x3)∈R3:0 ≤x1,x2,x3≤1}上的熱傳導(dǎo)問(wèn)題，如圖1 所示.邊界條件如下

這是文獻(xiàn)[11]采用的算例，文中沒(méi)有給出邊界量的計(jì)算結(jié)果，因而這里不便于比較.立方體的邊界劃分為54 個(gè)四邊形線性單元，每一正方形表面劃分9 個(gè)單元，如圖1 所示.圖2 給出了x3=1 面上的直線x1=11/18 上的9 個(gè)邊界節(jié)點(diǎn)處的溫度u和熱通量q1=?u/?x1的數(shù)值計(jì)算結(jié)果，可看出，數(shù)值解和精確解吻合的很好.此外，在立方體內(nèi)部選取400 個(gè)均勻分布在區(qū)域{0.15 ≤x1,x2≤0.85,x3=0.5}上的計(jì)算點(diǎn)，圖3(a)和圖3(b)分別展示了400 個(gè)計(jì)算點(diǎn)上的溫度u和通量q1=?u/?x1數(shù)值解的相對(duì)誤差曲面.

圖1 邊界網(wǎng)格及邊界計(jì)算點(diǎn)Fig.1 Boundary grid and boundary calculation points

圖2 圖1 中9 個(gè)邊界點(diǎn)處的溫度u 和通量q1Fig.2 Temperature and flux at 9 boundary points shown in Fig.1

圖3 場(chǎng)溫度和通量的相對(duì)誤差曲面Fig.3 Relative error surfaces for temperature and flux

算例2為了與文獻(xiàn)[11]提出的規(guī)則化邊界元法進(jìn)行比較，本算例采用文獻(xiàn)[11]中的算例.考慮單位球上的熱傳導(dǎo)，如圖4 所示.問(wèn)題具有Dirichlet 邊界條件：u(r,φ,θ)|Γ=，其中由下列問(wèn)題的解析解給出

在單位球內(nèi)部選擇400 個(gè)計(jì)算點(diǎn)，均勻分布在球面(根據(jù)θ,φ)上.圖5 描述了這400 個(gè)計(jì)算點(diǎn)上的場(chǎng)溫度u和場(chǎng)梯度?u/?x1數(shù)值解的L2誤差隨邊界單元數(shù)增加的變化情況，即收斂曲線.表1 列出了本文方法與規(guī)則化邊界積分方程法[11]在不同單元數(shù)下，求得的邊界位勢(shì)梯度q=?u/?x1和法向梯度q=?u/?n數(shù)值解的相對(duì)誤差REq1和REq，從對(duì)比中可看出，輔助邊值問(wèn)題法比規(guī)則化邊界積分方程的準(zhǔn)確性稍好，這主要是因?yàn)檩o助邊值問(wèn)題法計(jì)算系數(shù)矩陣的對(duì)角元可能更準(zhǔn)確.

圖4 單位球域Fig.4 Unit sphere

圖5 場(chǎng)溫度u 和通量q1=?u/?x1的收斂曲線Fig.5 Relative errors for the temperature and its derivative vs.the number of boundary element

表1 不同單元數(shù)下，REq 和REq1的數(shù)值結(jié)果Table 1 The numerical results of REq and REq1with different boundary elements

算例3在最后一個(gè)算例中，考慮一個(gè)復(fù)雜區(qū)域上的熱傳導(dǎo)問(wèn)題.計(jì)算區(qū)域是如圖6 所示的一個(gè)變形長(zhǎng)方體[31]，它是由長(zhǎng)方體{L×W×H=5.0 ×1.0×1.0}變形得到，其中L，W，H分別代表長(zhǎng)方體的長(zhǎng)、寬、高.變形長(zhǎng)方體的上、下面的兩條長(zhǎng)邊可分別表示為，左、右側(cè)面的寬度和高度不變.在變形長(zhǎng)方體邊界上施加混合邊界條件，其中變形長(zhǎng)方體下面通量已知，其余各面溫度已知，問(wèn)題的解析解為u(x1,x2,x3)=x1x2x3+10x1+10x2+10x3.

圖6 變形長(zhǎng)方體邊界網(wǎng)格Fig.6 The boundary meshes of the deformed cuboid

采用混合單元計(jì)算.變形長(zhǎng)方體的邊界被劃分為48 個(gè)單元，其中左、右側(cè)面各劃分為4 個(gè)線性單元(精確單元)，其余4 個(gè)面各劃分為10 個(gè)8 節(jié)點(diǎn)四邊形2 次單元，如圖5 所示.沿著變形長(zhǎng)方體的中心線選取10 個(gè)計(jì)算點(diǎn)，表2 列出了10 個(gè)計(jì)算點(diǎn)上的溫度u數(shù)值解以及精確解.此外，圖7 給出了變形長(zhǎng)方體的5 個(gè)表面(不包括下表面)上的邊界節(jié)點(diǎn)處的法向梯度?u/?n和梯度?u/?x1的相對(duì)誤差隨邊界單元數(shù)增大的變化情況，即收斂曲線.

表2 沿著長(zhǎng)方體中心線上的點(diǎn)的溫度u 的數(shù)值結(jié)果Table 2 Numerical results of the temperature of a point along the center line of a cuboid

圖7 邊界法向梯度?u/?n 和梯度?u/?x1的收斂曲線Fig.7 Relative errors for the boundary quantities and derivative vs.the number of boundary elements

4 結(jié)論

本文提出求解位勢(shì)梯度邊界積分方程的輔助邊值問(wèn)題法.構(gòu)造與邊值問(wèn)題具有相同解域的輔助邊值問(wèn)題，通過(guò)求解輔助邊值問(wèn)題，可準(zhǔn)確地獲得梯度邊界積分方程的離散系統(tǒng)矩陣，計(jì)算過(guò)程不涉及強(qiáng)奇異邊界積分的計(jì)算.所求得的系統(tǒng)矩陣可直接用來(lái)求解邊值問(wèn)題，不再需要重新計(jì)算系統(tǒng)矩陣.三個(gè)數(shù)值算例驗(yàn)證了方法的可行性、準(zhǔn)確性及收斂性.此外，從與規(guī)則化邊界元法的數(shù)值結(jié)果的對(duì)比可看出，本文方法的數(shù)值結(jié)果稍好一點(diǎn)，說(shuō)明采用輔助邊值問(wèn)題法計(jì)算梯度邊界積分方程的系統(tǒng)矩陣更準(zhǔn)確，特別是對(duì)角元的計(jì)算.

本文為梯度邊界積分方程的求解提供了新的思路.輔助邊值問(wèn)題方法具有理論簡(jiǎn)單，程序設(shè)計(jì)容易，精度高等優(yōu)點(diǎn)，而且容易拓展到彈性問(wèn)題、Stokes問(wèn)題、Helmholtz問(wèn)題等.