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

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

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

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

虛邊界元法是一種無奇異邊界元法，其最大優(yōu)點(diǎn)是能夠避免像邊界元法求解邊界物理量時(shí)遇到的奇異邊界積分和求解近邊界點(diǎn)處的物理量時(shí)遇到的幾乎奇異邊界積分.至今，虛邊界元法已廣泛應(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)的研究開辟了新的途徑，同時(shí)也拓展了虛邊界法的應(yīng)用領(lǐng)域.

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

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

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

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)分別為場點(diǎn)和源點(diǎn)，

令

則r(x,y)可表示為

r(x,y)=

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

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

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

(6)

(7)

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

(8)

(9)

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

(10)

(11)

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

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

GΦ=F

(12)

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

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

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

3 數(shù)值算例

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

(13)

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

算例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個(gè)單元，兩長邊各劃分為10個(gè)單元，共劃分為24個(gè)單元.在薄體矩形的兩短邊邊界上各均勻地選取20個(gè)計(jì)算點(diǎn)；在兩長邊邊界上各均勻選取40個(gè)計(jì)算點(diǎn)，虛、實(shí)邊界間的距離d分別取1、20、40，計(jì)算不同厚度h的板在短邊邊界上40個(gè)點(diǎn)處的溫度和長邊邊界上80個(gè)點(diǎn)處的通量.圖2和圖3描述了所得結(jié)果的平均相對誤差.從中可看出，當(dāng)板的厚度從10-1m到10-9m變化時(shí)，d=20和40時(shí)的數(shù)值結(jié)果具有很高的精度；d=1時(shí)，結(jié)果的精度稍差，但仍具有較高的精度.這表明虛、實(shí)邊界間的距離選取相當(dāng)?shù)貙挿?

圖2 邊界點(diǎn)溫度的平均相對誤差

圖3 邊界點(diǎn)通量的平均相對誤差

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

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

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

圖5 邊界點(diǎn)溫度和通量的收斂曲線

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

圖7 外邊界點(diǎn)溫度的平均相對誤差

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

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

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

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

圖10 邊界點(diǎn)溫度和通量的收斂曲線

4 結(jié)束語

本文提出求解二維各向異性位勢薄體問題的虛邊界元法，給出求解此類問題的新途徑，同時(shí)拓展了虛邊界元法的應(yīng)用領(lǐng)域.數(shù)值算例表明，虛邊界元法能夠非常有效地求解二維各向異性位勢薄體問題，而且虛、實(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] 張耀明，劉召顏，李功勝，等.各項(xiàng)異性位勢平面問題的規(guī)則化邊界元法[J].力學(xué)學(xué)報(bào), 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] 孫煥純，張立洲，許強(qiáng)，等.無奇異邊界元法[M].大連：大連理工大學(xué)出版社，1999.

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

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