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

李 聰 胡 斌 胡宗軍 牛忠榮

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

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

引言

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

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

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

1 二維位勢問題的常規邊界積分方程

考慮二維正交各向異性材料位勢問題,見圖1.內點y的位勢?(y)可由邊界位勢?(x)和法向位勢梯度q(x)的積分形式表達

式中,x(x1,x2)稱為場點,y(y1,y2)為源點.當源點y落在邊界Γ 上,可得二維正交各向異性位勢問題的邊界積分方程

圖1 正交各向異性材料位勢問題Fig.1 Potential problems of orthotropic materials

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

當源點y與積分單元Γe的節點重合,式(7) 的積分核是奇異的.當源點y臨近積分單元Γe,式(7)和式(8)的積分核呈現出幾乎奇異性.相較于奇異積分,高階單元幾乎奇異積分的計算更加復雜.為了盡可能減少處理高階單元奇異和幾乎奇異積分所帶來的計算量負擔,在復數域內建立高階單元奇異和幾乎奇異積分的計算列式.

2 復勢法計算高階單元奇異和近奇異積分

圖2 復平面坐標轉換Fig.2 Complex plane coordinate conversion

式中,Re[···]和Im[···]分別表示復變函數的實部和虛部.

2.1 線性單元

對于線性單元,如圖3 所示,z1,z2是單元首、末節點復數坐標.在單元局部坐標系oξ 下有

圖3 線單元坐標變換Fig.3 Coordinate conversion of a linear element

2.1.1 奇異積分

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

2.2 3 節點二次單元

對于3 節點二次單元,如圖4 所示,z1是單元首節點的復數坐標,z2是單元中節點的復數坐標,z3是單元末節點的復數坐標.在局部坐標系oξ 下有

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

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

2.2.1 情形Iza0

(1)奇異積分

當源點zs與積分單元首節點z1重合,使用ξ=2η ?1 變換可得

式(41)等號右邊第1 項采用常規Gauss 數值積分計算,第2 項采用對數型Gauss 數值積分計算.

當源點zs與積分單元中節點z2重合

式(42)等號右邊第1 項采用常規Gauss 數值積分計算,后2 項采用對數型Gauss 數值積分計算.

當源點zs與積分單元末節點z3重合,使用ξ=1 ?2η 變換可得

式(43)等號右邊第1 項采用常規Gauss 數值積分計算,第2 項采用對數型Gauss 數值積分計算.

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

當源點zs位于積分單元外

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

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

2.2.2 情形IIza=0

當za=0 時,如圖5 所示,式(32)變成

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

至此,線性單元和二次單元近場積分的計算已在復平面內完成.值得指出,復平面的使用極大地簡化了積分的計算公式,這為高階單元應用于快速多極邊界元法提供了便利.

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

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

圖6 “近場積分”和“遠場積分”Fig.6 “The near-field integrals”and“the far-field integrals”

3 算例

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

3.1 正交各向異性不規則閉域熱傳導問題

不規則閉域的外邊界為

情況1:當傳熱系數?1=?2時,各向同性問題.分別采用FM-BEM-C,FM-BEM-LA 和FM-BEM-QSA對圖7 結構進行計算,計算結果見表1.

圖7 不規則閉域熱傳導問題Fig.7 Heat conduction in an irregular domain

表1 不規則閉域A、B 點的溫度?Table 1 The temperatures of inner points A and B in an irregular domain

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

情況2:當傳熱系數?1=3,?2=2,正交各向異性問題.沿外邊界均勻劃分2880 個單元,沿內邊界均勻劃分720 個單元.采用FM-BEM-QSA 和FM-BEMLA 計算邊界C.D點附近的溫度? 和溫度梯度qx2,如表2 和表3 所示.

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

3.2 正交各向異性橢圓環閉域熱傳導問題

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

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

圖8 橢圓環閉域熱傳導問題Fig.8 Heat conduction in the elliptic ring domain

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

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

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

圖10 源點y 相對單元幾何說明Fig.10 Geometric relation between the source point y and element

3.3 隨機多孔方板

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

采用FM-BEM 對隨機多孔方板進行計算,通過調整隨機橢圓孔的數量,獲得采用FM-BEM-C,FMBEM-LA 和FM-BEM-QSA 消耗的CPU 時間與自由度(DOFs)之間的關系,如圖12 所示.

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

圖12 CPU 時間與自由度的關系Fig.12 Relationship between CPU time and degrees of freedom

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

4 結論

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

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

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

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