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

陳關(guān)忠， 周愛華， 張耀明

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

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

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

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

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

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

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

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

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

(1)

邊界條件為

(2)

(3)

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

其通量為

(4)

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

(5)

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

令

則r(x,y)可表示為

r(x,y)=

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

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

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

(6)

(7)

計算內(nèi)點位勢和位勢梯度的邊界積分方程表示為

(8)

(9)

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

(10)

(11)

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

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

GΦ=F

(12)

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

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

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

3 數(shù)值算例

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

(13)

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

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

圖1 各向異性薄體矩形區(qū)域的熱傳導(dǎo)

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

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

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

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

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

圖4 內(nèi)點溫度的相對誤差曲面

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

圖6 各向異性薄體圓環(huán)區(qū)域的熱傳導(dǎo)

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

圖8 內(nèi)邊界點通量的平均相對誤差

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

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

圖9 內(nèi)點溫度的相對誤差曲面

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

4 結(jié)束語

本文提出求解二維各向異性位勢薄體問題的虛邊界元法，給出求解此類問題的新途徑，同時拓展了虛邊界元法的應(yīng)用領(lǐng)域.數(shù)值算例表明，虛邊界元法能夠非常有效地求解二維各向異性位勢薄體問題，而且虛、實邊界距離的選取相當(dāng)?shù)貙挿?，即使結(jié)構(gòu)的厚度小到10-9m，依然可獲得高精度的數(shù)值解.

[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] 張耀明，劉召顏，李功勝，等.各項異性位勢平面問題的規(guī)則化邊界元法[J].力學(xué)學(xué)報, 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].大連：大連理工大學(xué)出版社，1999.

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

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