二維正交各向異性位勢(shì)問(wèn)題的高階單元快速多極邊界元法1)

李 聰 胡 斌 胡宗軍 牛忠榮

?(安徽建筑大學(xué)土木工程學(xué)院,合肥 230601)

?(合肥工業(yè)大學(xué)土木與水利工程學(xué)院力學(xué)系,合肥 230009)

引言

相較于有限元法(FEM),邊界元法[1](BEM) 只需在結(jié)構(gòu)外邊界劃分網(wǎng)格,因此BEM 可對(duì)復(fù)雜幾何模型進(jìn)行簡(jiǎn)單的網(wǎng)格劃分.然而BEM 生成的系數(shù)矩陣是非對(duì)稱滿陣,構(gòu)造矩陣所需的運(yùn)算量級(jí)為O(N2),如果采用高斯消去法求解,則所需運(yùn)算量級(jí)更是達(dá)到O(N3),因而導(dǎo)致BEM 難以計(jì)算大規(guī)模問(wèn)題.伴隨快速多極算法[2-3]發(fā)展起來(lái)的快速多極邊界元法(FMBEM)可很好地解決這個(gè)問(wèn)題,其所需的運(yùn)算量和存儲(chǔ)量較傳統(tǒng)BEM 大大降低,FM-BEM 成為近20 年邊界元領(lǐng)域的熱門課題之一,并被廣泛應(yīng)用于二維位勢(shì)和彈性問(wèn)題、三維聲場(chǎng)問(wèn)題等.

有關(guān)快速多極邊界元法的研究和應(yīng)用起源于20世紀(jì)90 年代初,主要尋求多極展開和局部展開概念在不同力學(xué)問(wèn)題的具體實(shí)現(xiàn)格式[4-6],而后被廣泛應(yīng)用于二維[7-14]、三維[15-22]問(wèn)題中.Pierce 等[10]采用二維實(shí)數(shù)域的Taylor 展開格式,實(shí)現(xiàn)了O(NlogN)量級(jí)的用于二維彈性力學(xué)問(wèn)題的多極邊界元法.Liu 等[11]將快速多極算法與對(duì)偶邊界元法相結(jié)合,模擬了二維線彈性多裂紋擴(kuò)展路徑.吳清華[12]將快速多極算法和高振蕩積分方法相結(jié)合,求解了一類帶高振蕩Hankel 核的邊界積分方程,求解結(jié)果和傳統(tǒng)方法吻合較好,計(jì)算效率也得到大幅度提升.Nishimura 等[19]采用分段常值形函數(shù)去離散裂紋問(wèn)題的超奇異邊界積分方程,然后耦合快速多極算法和GMRES 法求解三維裂紋問(wèn)題.趙麗濱等[20]利用Taylor 級(jí)數(shù)將三維彈性靜力學(xué)邊界積分方程的核函數(shù)進(jìn)行展開,引入快速多極算法應(yīng)用于薄體結(jié)構(gòu)問(wèn)題.Zhang 等[22]耦合了快速多極算法和混合邊界元節(jié)點(diǎn)法(Hybrid BNM),并將該方法成功應(yīng)用于三維位勢(shì)問(wèn)題的求解.Wang等[21]耦合快速多極邊界元法和重復(fù)相似子域法,構(gòu)造出含夾雜的復(fù)合材料結(jié)構(gòu)的預(yù)處理技術(shù),對(duì)三維復(fù)合材料進(jìn)行了大規(guī)模的數(shù)值模擬.

目前,快速多極邊界元法大多采用常值單元,其優(yōu)點(diǎn)是所有的近場(chǎng)積分,包括奇異積分和幾乎奇異積分[23-26]均可采用解析公式直接計(jì)算,但對(duì)于曲線邊界和復(fù)雜物理場(chǎng),常值單元的計(jì)算只能依賴于使用大量單元,且計(jì)算精度和計(jì)算效率往往不能令人滿意.使用高階單元能夠獲得更準(zhǔn)確的計(jì)算結(jié)果,但由于存在奇異積分和幾乎奇異積分的計(jì)算難題[27-28],因而鮮見(jiàn)高階單元快速多極邊界元法的相關(guān)報(bào)導(dǎo).本文在復(fù)平面內(nèi)計(jì)算高階單元的奇異積分和幾乎奇異積分,建立高階單元快速多極邊界元法,并將其應(yīng)用于二維正交各向異性位勢(shì)熱傳導(dǎo)問(wèn)題分析.

1 二維位勢(shì)問(wèn)題的常規(guī)邊界積分方程

考慮二維正交各向異性材料位勢(shì)問(wèn)題,見(jiàn)圖1.內(nèi)點(diǎn)y的位勢(shì)?(y)可由邊界位勢(shì)?(x)和法向位勢(shì)梯度q(x)的積分形式表達(dá)

式中,x(x1,x2)稱為場(chǎng)點(diǎn),y(y1,y2)為源點(diǎn).當(dāng)源點(diǎn)y落在邊界Γ 上,可得二維正交各向異性位勢(shì)問(wèn)題的邊界積分方程

圖1 正交各向異性材料位勢(shì)問(wèn)題Fig.1 Potential problems of orthotropic materials

其中,c(y)=α/(2π),α 為源點(diǎn)y處邊界的內(nèi)折角.位勢(shì)G(x,y)和位勢(shì)梯度F(x,y)的基本解為

當(dāng)源點(diǎn)y與積分單元Γe的節(jié)點(diǎn)重合,式(7) 的積分核是奇異的.當(dāng)源點(diǎn)y臨近積分單元Γe,式(7)和式(8)的積分核呈現(xiàn)出幾乎奇異性.相較于奇異積分,高階單元幾乎奇異積分的計(jì)算更加復(fù)雜.為了盡可能減少處理高階單元奇異和幾乎奇異積分所帶來(lái)的計(jì)算量負(fù)擔(dān),在復(fù)數(shù)域內(nèi)建立高階單元奇異和幾乎奇異積分的計(jì)算列式.

2 復(fù)勢(shì)法計(jì)算高階單元奇異和近奇異積分

圖2 復(fù)平面坐標(biāo)轉(zhuǎn)換Fig.2 Complex plane coordinate conversion

式中,Re[···]和Im[···]分別表示復(fù)變函數(shù)的實(shí)部和虛部.

2.1 線性單元

對(duì)于線性單元,如圖3 所示,z1,z2是單元首、末節(jié)點(diǎn)復(fù)數(shù)坐標(biāo).在單元局部坐標(biāo)系oξ 下有

圖3 線單元坐標(biāo)變換Fig.3 Coordinate conversion of a linear element

2.1.1 奇異積分

2.1.2 非奇異積分和幾乎奇異積分

2.2 3 節(jié)點(diǎn)二次單元

對(duì)于3 節(jié)點(diǎn)二次單元,如圖4 所示,z1是單元首節(jié)點(diǎn)的復(fù)數(shù)坐標(biāo),z2是單元中節(jié)點(diǎn)的復(fù)數(shù)坐標(biāo),z3是單元末節(jié)點(diǎn)的復(fù)數(shù)坐標(biāo).在局部坐標(biāo)系oξ 下有

注意到在使用二次單元模擬直線邊界時(shí),za可能等于0.因此,二次單元邊界積分的計(jì)算分為za=0 和za0 兩種情形.

圖4 二次單元坐標(biāo)變換Fig.4 Coordinate conversion of a three-node quadratic element

2.2.1 情形Iza0

(1)奇異積分

當(dāng)源點(diǎn)zs與積分單元首節(jié)點(diǎn)z1重合,使用ξ=2η ?1 變換可得

式(41)等號(hào)右邊第1 項(xiàng)采用常規(guī)Gauss 數(shù)值積分計(jì)算,第2 項(xiàng)采用對(duì)數(shù)型Gauss 數(shù)值積分計(jì)算.

當(dāng)源點(diǎn)zs與積分單元中節(jié)點(diǎn)z2重合

式(42)等號(hào)右邊第1 項(xiàng)采用常規(guī)Gauss 數(shù)值積分計(jì)算,后2 項(xiàng)采用對(duì)數(shù)型Gauss 數(shù)值積分計(jì)算.

當(dāng)源點(diǎn)zs與積分單元末節(jié)點(diǎn)z3重合,使用ξ=1 ?2η 變換可得

式(43)等號(hào)右邊第1 項(xiàng)采用常規(guī)Gauss 數(shù)值積分計(jì)算,第2 項(xiàng)采用對(duì)數(shù)型Gauss 數(shù)值積分計(jì)算.

(2)非奇異積分和幾乎奇異積分

當(dāng)源點(diǎn)zs位于積分單元外

二次單元坐標(biāo)變換的雅可比|J(ξ)|是根號(hào)下ξ 的二次多項(xiàng)式.根號(hào)的存在使得式(44)和式(45)難以處理.因此,將|J(ξ)|在ξ=0 處進(jìn)行泰勒展開

式(49)和式(50)等號(hào)右邊第1 項(xiàng)的奇異性已然去除,可采用常規(guī)的Gauss 數(shù)值積分計(jì)算.式(49)和式(50)等號(hào)右邊第2 項(xiàng)與式(46)和式(47)的積分項(xiàng)可歸納為以下幾種

2.2.2 情形IIza=0

當(dāng)za=0 時(shí),如圖5 所示,式(32)變成

坐標(biāo)變換的雅可比|J(ξ)| 和單元單位法向量n(ξ) 均為常量,此時(shí)單元是幾何線性的,奇異和幾乎奇異積分的計(jì)算可參考線性單元.

至此,線性單元和二次單元近場(chǎng)積分的計(jì)算已在復(fù)平面內(nèi)完成.值得指出,復(fù)平面的使用極大地簡(jiǎn)化了積分的計(jì)算公式,這為高階單元應(yīng)用于快速多極邊界元法提供了便利.

圖5 二次單元直邊界Fig.5 Three-node quadratic elements modelling the straight boundary

對(duì)于二維問(wèn)題,快速多極邊界元法利用四叉樹結(jié)構(gòu)將源點(diǎn)y對(duì)邊界單元的積分分成近場(chǎng)積分和遠(yuǎn)場(chǎng)積分.以圖6 為例,若源點(diǎn)位于cellC內(nèi),可將與cellC有公共角點(diǎn)的cell 稱之為“鄰居”.源點(diǎn)對(duì)于cellC及其“鄰居”內(nèi)的邊界單元的積分是近場(chǎng)積分,采用本文建立的算法處理.而源點(diǎn)對(duì)于剩余的邊界單元的積分是遠(yuǎn)場(chǎng)積分,采用快速多極展開式計(jì)算.有關(guān)快速多極展開式的計(jì)算過(guò)程這里不再闡述,詳見(jiàn)文獻(xiàn)[29].通過(guò)方程(2)求得所有未知的邊界位勢(shì)和位勢(shì)梯度后,域內(nèi)任意點(diǎn)的位勢(shì)和位勢(shì)梯度由式(1)和式(6)算得.

圖6 “近場(chǎng)積分”和“遠(yuǎn)場(chǎng)積分”Fig.6 “The near-field integrals”and“the far-field integrals”

3 算例

為驗(yàn)證本文建立的高階單元快速多極邊界元法的計(jì)算精度和計(jì)算效率,分別采用常值元快速多極邊界元法(FM-BEM-C)、線性元快速多極邊界元法(FM-BEM-LA)、二次元快速多極邊界元法(FM-BEMQSA)和二次元常規(guī)邊界元法(CBEM)計(jì)算了正交各向異性不規(guī)則閉域、圓環(huán)閉域熱傳導(dǎo)問(wèn)題和隨機(jī)多孔方板問(wèn)題.

3.1 正交各向異性不規(guī)則閉域熱傳導(dǎo)問(wèn)題

不規(guī)則閉域的外邊界為

情況1:當(dāng)傳熱系數(shù)?1=?2時(shí),各向同性問(wèn)題.分別采用FM-BEM-C,FM-BEM-LA 和FM-BEM-QSA對(duì)圖7 結(jié)構(gòu)進(jìn)行計(jì)算,計(jì)算結(jié)果見(jiàn)表1.

圖7 不規(guī)則閉域熱傳導(dǎo)問(wèn)題Fig.7 Heat conduction in an irregular domain

表1 不規(guī)則閉域A、B 點(diǎn)的溫度?Table 1 The temperatures of inner points A and B in an irregular domain

表1 顯示,若邊界單元總數(shù)n定值,與精確解相比,FM-BEM-QSA 的精度最高,FM-BEM-LA 次之,FM-BEM-C 最差.當(dāng)n=12 960 時(shí),B點(diǎn)FM-BEM-C的溫度? 計(jì)算值較解析解的相對(duì)誤差為0.77%.然而,在n=3240 時(shí),FM-BEM-LA 的計(jì)算值較解析解的相對(duì)誤差為0.15%;在n=1620 時(shí),FM-BEM-QSA 的計(jì)算值較解析解的相對(duì)誤差已小于0.12%.顯然,高階單元FM-BEM 的計(jì)算精度和計(jì)算效率較FM-BEM-C大幅度提升.

情況2:當(dāng)傳熱系數(shù)?1=3,?2=2,正交各向異性問(wèn)題.沿外邊界均勻劃分2880 個(gè)單元,沿內(nèi)邊界均勻劃分720 個(gè)單元.采用FM-BEM-QSA 和FM-BEMLA 計(jì)算邊界C.D點(diǎn)附近的溫度? 和溫度梯度qx2,如表2 和表3 所示.

表2 和表3 的結(jié)果顯示,FM-BEM-LA 和FMBEM-QSA 計(jì)算外邊界C.內(nèi)邊界D點(diǎn)附近的溫度?、溫度梯度值q與解析解吻合得很好,FM-BEM-QSA的計(jì)算結(jié)果更優(yōu).實(shí)際上,當(dāng)內(nèi)點(diǎn)十分接近外邊界C.內(nèi)邊界D點(diǎn)時(shí),近邊界點(diǎn)溫度的計(jì)算涉及高階單元的幾乎強(qiáng)奇異積分(1/r),溫度梯度的計(jì)算更是涉及幾乎超奇異積分(1/r2),計(jì)算難度較大,但本文第2 節(jié)建立的正則化算法能夠很好地處理這些積分.

3.2 正交各向異性橢圓環(huán)閉域熱傳導(dǎo)問(wèn)題

表2 正交各向異性不規(guī)則閉域C 點(diǎn)附近沿x2=0 mm 內(nèi)點(diǎn)的溫度? 和溫度梯度qx2Table 2 Temperature and temperature gradient of inner point near point C(x2=0 mm)

表3 正交各向異性不規(guī)則閉域D 點(diǎn)沿x1=10 附近的溫度? 和溫度梯度qx1Table 3 Temperature and temperature gradient of inner point near point D(x1=10 mm)

圖8 橢圓環(huán)閉域熱傳導(dǎo)問(wèn)題Fig.8 Heat conduction in the elliptic ring domain

定義厚長(zhǎng)比δ=(b2?b1)/b1,δ 越小,結(jié)構(gòu)越薄.使用FM-BEM 處理薄體問(wèn)題時(shí),近場(chǎng)積分的計(jì)算不僅僅涉及到奇異積分的計(jì)算,還涉及到幾乎奇異積分的計(jì)算.當(dāng)δ=0.02 時(shí),正交各向異性橢圓環(huán)全域的溫度云圖和y1方向的溫度梯度云圖如圖9 所示.表4 給出了取不同δ 的橢圓閉域結(jié)構(gòu)中A,B,C,D點(diǎn)的溫度? 和溫度梯度q的常規(guī)邊界元法(CBEM)和FM-BEM 計(jì)算值,這里CBEM(二次元) 采用的是二次單元離散邊界.與解析解相比,CBEM 和FM-BEM的計(jì)算結(jié)果在δ=1.0×10?1～1.0×10?6范圍內(nèi)均保持高精度.當(dāng)δ 減小至1.0×10?6時(shí),FM-BEM-LA 計(jì)算結(jié)果失效,CBEM(二次元) 和FM-BEM-QSA 的計(jì)算結(jié)果依舊精確,與解析解的相對(duì)誤差均小于0.1%.這是由于當(dāng)δ 減小至1.0×10?6,若邊界繼續(xù)采用線性元離散,源點(diǎn)y落在了域外,其計(jì)算結(jié)果失效,而若邊界采用二次元離散,源點(diǎn)y落在域內(nèi),其結(jié)果仍然精確,見(jiàn)圖10.因此在離散超薄結(jié)構(gòu)曲線型邊界時(shí),應(yīng)采用二次元.對(duì)于超薄體結(jié)構(gòu),FM-BEM-QSA 的計(jì)算精度與二次元CBEM 幾乎相同,在快速算法的幫助下,FM-BEM-QSA 能夠處理大規(guī)模問(wèn)題,具有更廣闊的應(yīng)用前景.

圖9 正交各向異性橢圓環(huán)閉域全域的位勢(shì)和位勢(shì)梯度云圖(δ=0.02)Fig.9 Distributions of potential and potential gradient of the orthotropic elliptic ring domain(δ=0.02)

表4 正交各向異性橢圓環(huán)閉域A,B,C,D 點(diǎn)的溫度? 和溫度梯度qTable 4 Temperature and temperature gradient of point A,B,C,D

圖10 源點(diǎn)y 相對(duì)單元幾何說(shuō)明Fig.10 Geometric relation between the source point y and element

3.3 隨機(jī)多孔方板

考慮一個(gè)邊長(zhǎng)200×200 含多個(gè)隨機(jī)孔的方板,如圖11 所示,m×m個(gè)隨機(jī)橢圓孔置于矩形方板中,其中m=100,200,300,400,500,600,800.傳熱系數(shù)?1=1,?2=2,每個(gè)橢圓孔邊界采用360 個(gè)節(jié)點(diǎn)離散,方板的每條直邊采用5000 個(gè)節(jié)點(diǎn)離散.矩形方板的左右兩直邊以及橢圓孔邊的位勢(shì)? 已知,矩形方板的上下兩直邊的位勢(shì)梯度q已知.該問(wèn)題的解析解為?=+3x1x2+5.所有的計(jì)算在一臺(tái)Intel Xeon(R) CPU E5-2670 @2.60GHz 內(nèi)存為16.0 GB 的電腦上完成.

采用FM-BEM 對(duì)隨機(jī)多孔方板進(jìn)行計(jì)算,通過(guò)調(diào)整隨機(jī)橢圓孔的數(shù)量,獲得采用FM-BEM-C,FMBEM-LA 和FM-BEM-QSA 消耗的CPU 時(shí)間與自由度(DOFs)之間的關(guān)系,如圖12 所示.

圖12 計(jì)算結(jié)果顯示FM-BEM-C,FM-BEM-LA 和FM-BEM-QSA 所需CPU 時(shí)間與自由度數(shù)量成線性關(guān)系,計(jì)算量均處于O(N)量級(jí).

圖12 CPU 時(shí)間與自由度的關(guān)系Fig.12 Relationship between CPU time and degrees of freedom

上述3 個(gè)算例的計(jì)算結(jié)果表明,在給定的計(jì)算精度下,相對(duì)于常值單元FM-BEM,高階單元FM-BEM使用節(jié)點(diǎn)數(shù)顯著減少,所需CPU 時(shí)間大幅度降低,因此高階單元FM-BEM 可更加有效求解大規(guī)模問(wèn)題.此外,高階單元奇異積分和幾乎奇異積分計(jì)算難題的解決,使得高階單元FM-BEM 能夠處理超薄體問(wèn)題,拓寬了高階單元FM-BEM 的適用范圍.

4 結(jié)論

通過(guò)解決高階單元奇異積分和幾乎奇異積分的計(jì)算難題,建立了適用于二維正交各向異性位勢(shì)問(wèn)題的高階單元快速多極邊界元(FM-BEM),其優(yōu)點(diǎn)如下:

(1) 高階單元FM-BEM 的計(jì)算精度和計(jì)算效率較常值單元FM-BEM 的精度高.

(2)高階單元奇異積分和幾乎奇異積分難題的解決,使得高階單元FM-BEM 能夠應(yīng)用于超薄體結(jié)構(gòu),具有廣泛的應(yīng)用前景.

(3) 高階單元FM-BEM 計(jì)算時(shí)間與自由度數(shù)量成線性關(guān)系,其計(jì)算效率仍處于O(N) 量級(jí),適用于求解大規(guī)模問(wèn)題.