董榮榮 張 超 張耀明

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

引言

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

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

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

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

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

1 問題描述

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

具有邊界條件

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

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

2 位勢梯度邊界積分方程的輔助邊值問題法

2.1 邊界積分方程

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

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

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

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

這里ξ1,ξ2是無因次坐標，且-1 ≤ξ1,ξ2≤1.因此單元上的點x可以近似地表示為

這里xk是第k個插值節(jié)點的坐標.

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

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

2.3 精確單元

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

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

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

這里θi,φi(i=1,2)是精確單元的端點坐標.

2.4 邊界積分方程的離散化

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

2.5 輔助邊值問題法

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

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

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

3 數(shù)值算例

考慮3 個數(shù)值算例，來驗證方法的有效性.數(shù)值實驗的重點在于考察輔助邊值問題法計算邊界通量?u/?x1的能力及準確性.為了估計數(shù)值誤差，采用如下L2范數(shù)

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

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

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

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

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

圖3 場溫度和通量的相對誤差曲面Fig.3 Relative error surfaces for temperature and flux

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

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

圖4 單位球域Fig.4 Unit sphere

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

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

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

表2 沿著長方體中心線上的點的溫度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é)論

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

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