二維各向異性位勢薄體問題的虛邊界元法

陳關忠， 周愛華， 張耀明

(山東理工大學 理學院， 山東 淄博 255091)

隨著科學技術的進步與發展，各向異性材料由于在結構和性能方面具有眾多的優越性而在工程中應用日趨廣泛.特別是各向異性薄膜材料、薄體結構，如機械構件涂層、復合材料疊層構件、薄殼類構件、靈敏構件、以及纖維增強材料的粘結層等已廣泛應用于工程中.但受其厚度尺寸的限制，結構物理量的數值分析一直是工程中的難點[1].采用有限元分析時，為了避免出現畸形單元，需按照薄體的厚度劃分單元.可是，這樣做必然會導致出現成千上萬個單元，計算工作量劇增，對此小型計算機甚至無法完成.

邊界元法能夠有效地求解各向同性薄體結構問題[2-3]，但是奇異基本解的使用導致邊界積分方程中不可避免地出現奇異邊界積分和幾乎奇異邊界積分，這些積分的處理和計算具有相當的難度且耗時.各向異性薄體結構問題的邊界元法是一個更困難的問題，至今取得的成果很少[4].

虛邊界元法是一種無奇異邊界元法，其最大優點是能夠避免像邊界元法求解邊界物理量時遇到的奇異邊界積分和求解近邊界點處的物理量時遇到的幾乎奇異邊界積分.至今，虛邊界元法已廣泛應用于各種常規結構[5]及各向同性薄體結構[6-7]，但對于各向異性狹窄結構的研究尚未發現.本文拓展了虛邊界法的應用范圍，將其應用于二維位勢各向異性狹窄結構，為各向異性薄體結構的研究開辟了新的途徑，同時也拓展了虛邊界法的應用領域.

虛邊界元法的求解精度受虛實邊界間的距離影響.本文選擇近、中、遠三種不同的虛、實邊界間的距離，對不同厚度的各向異性薄體問題進行了計算，均獲得了令人滿意的數值結果. 表明虛、實邊界間的距離選取非常地寬泛，且方法簡單、效率高，即使結構的寬度狹窄到納米級(10-9m)，依然可獲得很高精度的數值解.

1 二維各向異性位勢薄體問題的虛邊界元法

本文假定Ω是R2中的一個有界區域，Γ=?Ω是邊界.n=(n1,n2)是區域Ω的邊界Γ在x點處的單位外法向量.

1.1 二維各向異性位勢薄體問題的控制方程

二維各向異性位勢薄體問題的控制方程為

(1)

邊界條件為

(2)

(3)

式中u為勢函數，kij(i,j=1,2)為各向異性材料特性系數，n為邊界外法向量，Γ1,Γ2分別是u和?u/?n已知的邊界.

其通量為

(4)

二維各向異性位勢薄體問題控制方程的基本解為

(5)

其中，x=(x1,x2),y=(y1,y2)分別為場點和源點，

令

則r(x,y)可表示為

r(x,y)=

1.2 二維各向異性位勢薄體問題虛邊界積分方程

設想假定Ω被嵌入一個無限的區域中，在Γ外的無限區域中有一延拓邊界Γ′(這里稱為虛邊界)，沿著邊界Γ′有一個待定的虛擬密度函數Φ(y)，令此虛擬密度函數對真實邊界所產生的位勢或法向梯度與邊界條件一致，從而達到求解待定密度函數的目的.以上稱之為虛邊界元法.

各向異性位勢薄體問題中的虛邊界積分方程為

(6)

(7)

計算內點位勢和位勢梯度的邊界積分方程表示為

(8)

(9)

2 虛邊界積分方程的離散化方法

(10)

(11)

內點位勢和位勢梯度的離散公式可如法炮制，這里不再贅述.

不失一般性，假定邊界為Dirichlet邊界條件，則式(10)可以轉化為如下矩陣形式

GΦ=F

(12)

這里G=[Gij]n×n，Φ=[Φ1,Φ2,…，Φn]T,F=[u(x1),u(x2),…，u(xn)]T,其中

由于此方法的積分點和場點分別位于虛邊界和實邊界上，因此無需考慮奇異積分，Gij可直接用標準高斯積分進行計算.本文算例均采用八點高斯積分公式計算.

由方程組(12)可求解密度函數Φ(y)，然后代入式(7)-(9)，進而依次求得邊界的位勢法向梯度、內點的位勢和位勢梯度.

3 數值算例

我們考慮一個各向異性薄體矩形和一個各向異性薄體圓環的數值算例，來驗證本文方法的有效性，所有算例采用常單元等額配點法.為了表明方法數值解的準確性，定義平均相對誤差如下：

(13)

虛邊界元法的求解精度受虛實邊界間的距離影響.以下算例中，我們選擇近、中、遠三種不同的虛實邊界間的距離，對不同厚度的各向異性薄體問題進行了計算，都獲得了令人滿意的數值結果，表明虛實邊界間的有效距離選取范圍非常地寬泛.

算例1各向異性薄體矩形區域的熱傳導.如圖1所示，薄體矩形區域的長度l=1m，厚度為h從10-1m到10-9m之間變化，各向異性材料的系數為k11=2,k12=1,k22=3，邊界條件為在短邊上已知通量q，在長邊上已知溫度u.本問題的精確解為u(x,y)=x2-y2+xy.

圖1 各向異性薄體矩形區域的熱傳導

如圖1所示，取虛邊界與結構外表面的幾何形狀相似，將兩短邊各劃分為2個單元，兩長邊各劃分為10個單元，共劃分為24個單元.在薄體矩形的兩短邊邊界上各均勻地選取20個計算點；在兩長邊邊界上各均勻選取40個計算點，虛、實邊界間的距離d分別取1、20、40，計算不同厚度h的板在短邊邊界上40個點處的溫度和長邊邊界上80個點處的通量.圖2和圖3描述了所得結果的平均相對誤差.從中可看出，當板的厚度從10-1m到10-9m變化時，d=20和40時的數值結果具有很高的精度；d=1時，結果的精度稍差，但仍具有較高的精度.這表明虛、實邊界間的距離選取相當地寬泛.

圖2 邊界點溫度的平均相對誤差

圖3 邊界點通量的平均相對誤差

在薄體矩形區域內均勻地選取100個計算點，當h=10-9m和d=10時，計算這些點處的溫度，圖4給出了計算結果的誤差曲面. 可看出，數值結果的相對誤差非常地小.這表明該方法能夠有效地求解厚度小到納米級的各向異性薄體問題.

當h=10-9m和d=20時，隨著邊界單元數的增加，計算短邊邊界上40個點處的溫度和長邊邊界上80個點處的通量，其數值結果的平均相對誤差變化列于圖5.計算時，長、短邊劃分單元數的比例為5∶1.可看出，隨著單元數的增加，相對誤差迅速地減小，說明該方法具有良好的收斂性.

圖4 內點溫度的相對誤差曲面

圖5 邊界點溫度和通量的收斂曲線

圖6 各向異性薄體圓環區域的熱傳導

圖7 外邊界點溫度的平均相對誤差

圖8 內邊界點通量的平均相對誤差

在薄體圓環區域內均勻地選取100個計算點，當h=10-9m和d1=0.6,d2=20時，計算這些點處的溫度，圖9給出了計算結果的相對誤差曲面. 可看出，數值結果的相對誤差非常地小.這表明該方法能夠準確高效地求解厚度小到納米級的各向異性薄體問題.

當h=10-9m和d1=0.6,d2=20時，隨著邊界單元數的增加，計算外邊界上100個計算點處的溫度和內邊界上100個計算點處的通量，圖10給出了數值結果的收斂曲線.計算時，內外邊界劃分單元數的比例為1∶1.可看出，隨著單元數的增加，相對誤差迅速地減小，說明方法仍然具有較好的收斂性.

圖9 內點溫度的相對誤差曲面

圖10 邊界點溫度和通量的收斂曲線

4 結束語

本文提出求解二維各向異性位勢薄體問題的虛邊界元法，給出求解此類問題的新途徑，同時拓展了虛邊界元法的應用領域.數值算例表明，虛邊界元法能夠非常有效地求解二維各向異性位勢薄體問題，而且虛、實邊界距離的選取相當地寬泛，即使結構的厚度小到10-9m，依然可獲得高精度的數值解.

[1] Luo J F,Liu Y J,Berger E J.Analysis of two-dimensional thin structures(from micro-to nano-scales) using the boundary element method[J].Computational Mechanics,1998,22:404-412.

[2] Zhang Y M, Gu Y, Chen J T.Analysis of 2D thin walled structures in BEM with high-order geometry elements using exact integration[J]. Computer Modeling in Engineering and Sciences, 2009, 50(1):1-20.

[3] 張耀明，劉召顏，李功勝，等.各項異性位勢平面問題的規則化邊界元法[J].力學學報, 2011，43(4):785-788.

[4] Zhou H L, Niu Z R, Cheng C Z，etal. Analytical integral algorithm applied to boundary layer effect and thin body effect in BEM for anisotropic potential problems[J]. Computers and Structures,2008,86:1656-1671.

[5] 孫煥純，張立洲，許強，等.無奇異邊界元法[M].大連：大連理工大學出版社，1999.

[6]袁飛，屈文鎮，張耀明.二維薄體結構位勢問題的虛邊界元法[J].山東理工大學學報:自然科學版，2012,26(3)：18-21.

[7]屈文鎮，袁飛，張耀明.二維彈性薄體問題的虛邊界元法[J].山東理工大學學報：自然科學版，2012,26(3):64-67.