邊界元法的若干進(jìn)展和它在固體力學(xué)中的應(yīng)用

1、邊界元法的若干進(jìn)展和它在固體力學(xué)中的應(yīng)用 清華大學(xué)工程力學(xué)系清華大學(xué)工程力學(xué)系 姚振漢姚振漢引 言n彈性力學(xué)的三種提法微分提法變分提法積分提法偏微分方程邊值問題泛函極值問題邊界積分方程問題求解析解差分法求數(shù)值解李茲法求近似解有限元法有限元法求數(shù)值解邊界元法邊界元法求數(shù)值解歸結(jié)為求解線性代數(shù)方程組常化為常微分方程三種提法是完全等價(jià)的邊界元發(fā)展歷史回顧n邊界積分方程邊界元法有限元法（1955，56）之后發(fā)展起來的一種精確高效的工程數(shù)值分析方法在固體力學(xué)領(lǐng)域有限元法最重要的補(bǔ)充n邊界元法間接法位勢問題（Smith & Pierce, 1958)彈性力學(xué)（Massonet，1965）n邊界元法

2、直接法位勢問題（Jaswon，1963）彈性力學(xué)（Rizzo，1967）邊界元發(fā)展歷史回顧n19942003被SCI收錄的論文與邊界元法有關(guān)的有3904篇與有限元法有關(guān)的為16823篇n19902002被EI收錄的論文與邊界元法有關(guān)的有19968篇與有限元法有關(guān)的為75184篇與斷裂力學(xué)有關(guān)的為23647篇n在工程應(yīng)用方面在應(yīng)用最多的部門也從未超過有限元法的十分之一研究組邊界元研究歷史回顧n我們研究組邊界元法研究開始于1979年u基于彈性力學(xué)問題Rizzo型邊界積分方程邊界元法研究了彈性應(yīng)力集中問題和薄板彎曲問題u研究了邊界元有限元耦合方法u研究了邊界元法在形狀優(yōu)化缺陷識別等逆問題中的應(yīng)用uG

3、alerkin對稱邊界元法用于結(jié)構(gòu)極限與安定分析等問題u研究了精確高效的計(jì)算方法，提出了邊界元法誤差的一種直接估計(jì)。u對于彈性接觸問題，提出了單元與單元間協(xié)調(diào)的接觸方案，研究了二維、三維移動、滾動接觸。研究組邊界元研究近年工作u2000前后針對復(fù)合材料，對于含隨機(jī)分布大量夾雜的二維彈性固體提出了一種重復(fù)相似子域邊界元法，計(jì)算了100多個夾雜、近萬自由度的問題。u研究了在微機(jī)機(jī)群上的并行算法。2000年在由8臺微機(jī)組成的機(jī)群上最大計(jì)算規(guī)模45,000自由度。n近年來多極快速算法在邊界元法中的應(yīng)用給邊界元法解決復(fù)雜工程與科學(xué)問題展示了廣闊的前景。u用于含隨機(jī)分布夾雜二維、三維彈性體數(shù)值模擬，一臺微

4、機(jī)可計(jì)算數(shù)十萬自由度的問題，在微機(jī)機(jī)群并行系統(tǒng)最大的二維算例有800萬自由度。u計(jì)算了含16384條隨機(jī)分布裂紋的二維無限彈性體，1,572,864自由度，研究了相應(yīng)的裂紋擴(kuò)展問題。鑒于u 邊界元法始終是計(jì)算力學(xué)中一個重要的研究領(lǐng)域u 快速多極邊界元法是近年受到特別關(guān)注的一個研究方向u 本人已經(jīng)從事邊界元法研究25年，幾乎沒有間斷u 近年研究組研究重點(diǎn)放在快速多極邊界元法的研究方面因此今天就以此為題，希望能對各位同學(xué)的學(xué)習(xí)有一定的幫助。彈性力學(xué)邊界積分方程n彈性靜力學(xué)的邊界積分方程其中SSS(p)(p)=(p; Q)(Q)dV(Q)+(p; q) (q)dS(q) (p; q)(q)dS(q)

5、ijjijjijjVSijjSCuufuttuS1(P; Q)=(34 ),16(1)ijijijur rGr S21(P; q) =8 (1) (12 )3 . ,(12 )( ,)ijijijijjitrrr rr nrnn 對于同一彈性體同一彈性體的兩種變形狀態(tài)兩種變形狀態(tài)ftii( )( ),11 )1 ()1 ()1 ( , ,ijijiuftii( )( ),22 )2()2()2( , ,ijijiuViSiVufSutd+d)1(2)i)1(2)iVijijSijijVuSund,d)1(2)1(2)VijijVjiijVuVud,d,)1(2)1(2)VjiijVud,)1

6、(2)VijijVd)1 (2)VijklijklVEd)1 (2)VijijVd)2(1)ViSiVufSutd+d)2(1)i)2(1)i通常一個狀態(tài)是待求狀態(tài)通常一個狀態(tài)是待求狀態(tài)另一狀態(tài)是已知的輔助狀態(tài)另一狀態(tài)是已知的輔助狀態(tài)彈性力學(xué)的邊界積分方程可以由Betti定理定理出發(fā)導(dǎo)出以無限彈性空間中一點(diǎn)作用單位集中力的Kelvin問題的解為輔助解0) ,() ;(,) ;(,)(SSQPQPGuQPuGijkkijkikjjiijijrrrGQPu,)43()1 (161=) ;(S由Betti定理可得Somigliana等式（無體積力情況）dS(q)(q)q) P;(dS(q)(q)q)

7、 P;(=P)(SSSjijSjijiuttuu再將集中力作用點(diǎn)（源點(diǎn)）再將集中力作用點(diǎn)（源點(diǎn)）P 趨趨于邊界點(diǎn)于邊界點(diǎn) p，考慮到源點(diǎn)是積分，考慮到源點(diǎn)是積分的奇異點(diǎn)，作適當(dāng)處理，即可得的奇異點(diǎn)，作適當(dāng)處理，即可得到前面給出的到前面給出的邊界積分方程邊界積分方程。彈性力學(xué)邊界元法n將邊界分成邊界單元在每個單元上將邊界變量插值離散可得1eNjj ()()11()()11( )( )( ) ( )( )( )( )( ) ( ) ( )( ) ( )mmmmallallllmmmmallalllluNuuNutNttNt()11()11( )( , ( )( )d ( )() =( , ( )(

8、 )d ( )()ejeuNmSmljljlNmSmljljlCp utp qNuqup qNtq 將弱解代入方程得到誤差對于加權(quán)余量法的配點(diǎn)格式，權(quán)函數(shù)采用Dirac-delta函數(shù)，要求即得()11()11( )( )( , ( )( )d ( )() n( , ( )( )d ( )()ejeuNmSmljljlNmSmljljlpCp utp qNuqup qNtq ( ) ( , )d (p)=0npp q()0, =1, 2, , nqnN彈性力學(xué)邊界元法將方程寫成矩陣形式，得其中H，G矩陣元素由核函數(shù)與形函數(shù)在單元上的積分求得。將邊界變量列矢量 U，T 按未知量與給定量重新排列，

9、可得邊界元法的求解代數(shù)方程組即 HUGT AXBC AXC彈性力學(xué)邊界元法n核函數(shù)與形函數(shù)乘積在單元上的積分矩陣元素都是核函數(shù)與形函數(shù)乘積在單元上的積分，矩陣是滿陣。主要計(jì)算量就是計(jì)算這些積分，以及求解滿陣代數(shù)方程組。對于規(guī)模不太大的問題，計(jì)算積分的工作量是主要計(jì)算量。非奇異積分采用等精度高斯積分格式求積，高斯點(diǎn)數(shù)由精度要求對不同情況自動確定。當(dāng)核函數(shù)的源點(diǎn)落入積分單元時出現(xiàn)奇異積分，包括弱奇異積分和柯西主值積分。彈性力學(xué)邊界元法邊界元法的優(yōu)缺點(diǎn)n邊界元法與有限元法及其它數(shù)值方法相比較的優(yōu)缺點(diǎn)u優(yōu)點(diǎn)：高精度（由于采用了解析基本解）降低了維數(shù)，便于模擬復(fù)雜邊界形狀對于高梯度、甚至有奇異性的問題，

10、不僅有較高的精度，而且在同等精度條件下有較高效率適合于處理無限域、半無限域問題適合于處理彈性接觸等邊界條件非線性問題u缺點(diǎn)：適用范圍遠(yuǎn)沒有有限元法廣泛解題規(guī)模受限制（方程系數(shù)矩陣為滿陣）對于域內(nèi)方程非線性問題優(yōu)勢減弱或喪失邊界元法的特點(diǎn)u早期有的文章通過簡單算例高估了邊界元法的計(jì)算精度，聲稱能達(dá)到0.01%，甚至更高。這里要分兩類問題，一類是沒有離散誤差的問題，稱為簡單問題，離散插值的邊界變量能精確滿足給定邊界條件，此時能達(dá)到很高精度，例如106，甚至108。對于一般問題，有離散誤差，只有合理劃分足夠多的邊界單元，才能達(dá)到要求的計(jì)算精度。u薄板梁純彎曲問題，只要采用二次單元，就是一個簡單問題。

11、通常認(rèn)為細(xì)長薄板梁要劃分較多單元才能達(dá)到滿意精度，我們只用4個二次元，就達(dá)104精度。u有限元法中單元邊長比不能太大，相鄰單元尺寸也不能相差太大。邊界元法則并不受此限制，上述100:1純彎薄板梁相鄰單元長度比為100，只要保證積分精度等運(yùn)算精度，還是可以得到高精度的結(jié)果。u有限元法中高斯積分通常只用12，或22個高斯積分點(diǎn)。邊界元法常采用等精度高斯積分，根據(jù)給定積分精度要求來確定用多少高斯點(diǎn)。對于上述細(xì)長薄板梁的算例，源點(diǎn)在邊長為1的短邊中點(diǎn)、在長邊二次單元上的積分，需要5060個高斯點(diǎn)進(jìn)行高斯積分。邊界元法的特點(diǎn)u邊界元法在解出未知的邊界變量之后，通過域內(nèi)變量的邊界積分公式可以求得域內(nèi)任意點(diǎn)

12、的位移、應(yīng)變和應(yīng)力。這樣得到的域內(nèi)位移場精確滿足Navier方程，但是當(dāng)它趨于邊界時的極限一般說來和邊界給定量或邊界元法解出的邊界變量并不一致。這種差別是和離散誤差密切相關(guān)的。對于沒有離散誤差的簡單問題，域內(nèi)解趨于邊界的極限和邊界變量是一致的，只有很小的運(yùn)算誤差。這也導(dǎo)致在有離散誤差的情況下，在近邊界區(qū)求得的內(nèi)點(diǎn)變量有較大的誤差，即通常稱為“邊界效應(yīng)”。邊界元法的特點(diǎn)邊界元法模擬光彈實(shí)驗(yàn)u為了便于將得到的結(jié)果和工程界熟悉的光彈性實(shí)驗(yàn)結(jié)果比較，域內(nèi)應(yīng)力分布也畫出了等色線和等傾線圖。而且我們一再表明，在數(shù)值模擬已經(jīng)很成熟的今天，進(jìn)行光彈性實(shí)驗(yàn)這樣的物理模擬已經(jīng)沒有多大意義。因?yàn)樗膊皇菍?shí)測，光彈性

13、模型和實(shí)際問題相比引進(jìn)了許多誤差，得到的結(jié)果并不比數(shù)值模擬結(jié)果更接近真實(shí)情況；而且數(shù)值模擬比光彈性模擬更經(jīng)濟(jì)便捷。下圖是帶孔板單向拉伸的等色線圖（條紋等級反映主應(yīng)力差）。邊界元法模擬光彈實(shí)驗(yàn)邊界元法并行計(jì)算的工程應(yīng)用o邊界元法結(jié)合并行計(jì)算，在8臺微機(jī)組成的機(jī)群并行系統(tǒng)上計(jì)算了一些規(guī)模稍大的工程實(shí)際問題。例如：u石油鉆桿的偏心鉆挺，49818自由度，8臺微機(jī)，計(jì)算43小時41分鐘。u渤海石油公司“濱海109”海洋鋪管船的錨機(jī)滾筒軸，24444自由度，7臺微機(jī)，計(jì)算436分鐘。這是2000年博士生尹欣得到的結(jié)果，當(dāng)時文獻(xiàn)中最大計(jì)算規(guī)模10萬自由度。子域結(jié)點(diǎn)數(shù)單元數(shù)13358111825331211

14、04總體166066638子域子域 1的網(wǎng)格劃分的網(wǎng)格劃分子域子域 25的網(wǎng)格劃分的網(wǎng)格劃分偏心鉆挺模型（5萬自由度）有限元模型（16000結(jié)點(diǎn)，13000千單元）鉆挺頂部Mises應(yīng)力圖（邊界元法，單位：MPa）鉆挺頂部Mises應(yīng)力圖（邊界元法，單位：MPa）鉆挺頂部Mises應(yīng)力圖（有限元法，單位：MPa）鉆挺頂部Mises應(yīng)力圖（有限元法，單位：MPa）滾筒軸邊界元模型：5個子域，8184結(jié)點(diǎn)，1721單元滾筒軸有限元模型：13000結(jié)點(diǎn)，64000單元邊界元解：左端局部Mises應(yīng)力（單位：MPa）有限元解：左端局部Mises應(yīng)力（單位：MPa）邊界元解：右端局部Mises應(yīng)力（單位

15、：MPa）有限元解：右端局部Mises應(yīng)力（單位：MPa）在復(fù)合材料模擬中的應(yīng)用n針對復(fù)合材料數(shù)值模擬，研究組對于含隨機(jī)分布孔洞或夾雜的二維彈性體進(jìn)行了深入的研究。提出了重復(fù)相似子域邊界元法，提高了計(jì)算效率。分析計(jì)算了含100個圓孔、圓形夾雜、橢圓形夾雜、多種形狀夾雜的二維彈性體，得到了等效材料常數(shù)和體積比、夾雜與基底材料常數(shù)比等的關(guān)系?？紤]了夾雜和基底之間有界面層的情況。對此類問題，引進(jìn)快速算法之后，使計(jì)算規(guī)模大大提高，初步工作已經(jīng)從100多個夾雜增加到了1600個夾雜。在復(fù)合材料模擬中的應(yīng)用重復(fù)相似子域邊界元法在復(fù)合材料模擬中的應(yīng)用含100個隨機(jī)分布圓形夾雜的二維彈性體在復(fù)合材料模擬中的應(yīng)

16、用含100個隨機(jī)分布圓形夾雜的二維彈性體變形圖在復(fù)合材料模擬中的應(yīng)用含100個隨機(jī)分布圓形夾雜的二維彈性體應(yīng)力分布在復(fù)合材料模擬中的應(yīng)用含100個隨機(jī)分布橢圓形夾雜的二維彈性體在復(fù)合材料模擬中的應(yīng)用含100個隨機(jī)分布橢圓形夾雜二維彈性體的變形圖在復(fù)合材料模擬中的應(yīng)用含100個隨機(jī)分布橢圓形夾雜二維彈性體的應(yīng)力分布在復(fù)合材料模擬中的應(yīng)用含100個隨機(jī)分布不同形狀夾雜的二維彈性體在復(fù)合材料模擬中的應(yīng)用含100個隨機(jī)分布不同形狀夾雜二維彈性體的變形圖在復(fù)合材料模擬中的應(yīng)用含100個隨機(jī)分布不同形狀夾雜二維彈性體的應(yīng)力分布在復(fù)合材料模擬中的應(yīng)用含100個隨機(jī)分布圓形夾雜（有界面層）的二維彈性體0246

17、81012141618200.50.91.31.72.12.52.93.3inclusion-matrix modulus ratioEE0dilutie solutionM-T, IDD solutionsCCM, GSCM solutionsBEM soutions Volume fraction c=0.5在復(fù)合材料模擬中的應(yīng)用含圓形夾雜二維彈性體等效彈性模量和模量比的關(guān)系在復(fù)合材料模擬中的應(yīng)用0246810121416182000.511.522.533.5inclusion-matrix modulus ratioKK0dilution solutionCCM, GSCM, M-T,

18、 IDD solutionsBEM solution Volume fraction c=0.5含圓形夾雜二維彈性體等效體積模量和模量比的關(guān)系在復(fù)合材料模擬中的應(yīng)用0246810121416182000.40.81.21.622.42.83.2inclusion-matrix modulus ratioGG0dilute solutionM-T, IDD solutionsCCM, GSCM solutionsBEM solutions Volume fraction c=0.5含圓形夾雜二維彈性體等效剪切模量和模量比的關(guān)系在復(fù)合材料模擬中的應(yīng)用024681012141618200.20.50

19、.81.11.41.722.3inclusion-matrix modulus ratioEE0dilute solutionM-T, IDD solutionsCCM, GSCM solutionsBEM solution Volume fraction c=0.4含圓形夾雜二維彈性體等效彈性模量和模量比的關(guān)系快速多極算法n快速多極算法 (fast multipole method, FMM)最初針對問題最初針對問題：多粒子相互作用的有勢場 存儲和計(jì)算復(fù)雜度：O(N2)uBarnes, Hut, (1986) 樹結(jié)構(gòu)算法 存儲和計(jì)算復(fù)雜度：O(NlogN)uGreengard, Rokhli

20、n, (1986) 快速多極算法 存儲和計(jì)算復(fù)雜度：O(N)uHrycak T, Rokhlin V, (1998) 新型快速多極算法 存儲和計(jì)算復(fù)雜度：O(N) 快速多極邊界元法n常規(guī)邊界元法由于方程組的系數(shù)矩陣為非對稱滿陣，對于N自由度系統(tǒng)，其運(yùn)算量為N 3量級，對于存儲的要求為N 2量級。邊界元快速算法通過快速多極展開，使運(yùn)算量和存儲量均減少到NlogN 量級。最近引入了局部展開的思想，進(jìn)一步使運(yùn)算量和存儲量要求減少到了N 量級。研究組初步的工作已經(jīng)把解題規(guī)模明顯擴(kuò)大，由原來8臺微機(jī)并行最多計(jì)算49,818自由度問題，到用邊界元快速算法一臺微機(jī)計(jì)算544,000自由度和1,572,864

21、自由度。 Sandia BEMSandia BEM軟件軟件(1994)(1994)Coyote AutoMEMS Coyote AutoMEMS 軟軟 件件評評 論論計(jì)算機(jī)CPU算法計(jì)算能力計(jì)算速度大規(guī)模并行機(jī)1900個Intel i860傳統(tǒng)邊界元法50,000個常值面元每小時140000常值元微機(jī)1個Intel 奔騰III多極快速邊界元法20,000,000個常值面元每小時1,200,000常值元 便宜了2000倍關(guān)鍵是改進(jìn)了算法規(guī)模增大400倍提高速度9倍邊界元快速算法和常規(guī)邊界元法的比較快速多極邊界元法基本思想n快速多極算法的基本思想：如果每個人要自己把信送到每個收信人，是非常費(fèi)時的，

22、當(dāng)收信人的范圍很大時是不可能的。但是現(xiàn)代的郵政系統(tǒng)，使每個人可以方便地給世界各地的收信人發(fā)信。傳統(tǒng)邊界元法形成滿陣方程組的系數(shù)要一個個地獨(dú)立計(jì)算，就像每個人親自送信一樣?？焖俣鄻O算法就像現(xiàn)代郵政系統(tǒng)一樣大大提高了效率。樹型存儲結(jié)構(gòu)n將求解區(qū)域與邊界有關(guān)部分劃分為四叉樹（二維問題）或八叉樹（三維問題）近場與遠(yuǎn)場對源點(diǎn) x，近場按常規(guī)邊界元法計(jì)算方程系數(shù)，遠(yuǎn)場用快速多極展開，不再計(jì)算一個個系數(shù)?？焖俣鄻O算法基本步驟n對于遠(yuǎn)場的快速多極算法有四個基本步驟：1.核函數(shù)的多極展開 多極展開系數(shù)（多極矩）g (y, k) 對所有核函數(shù)源點(diǎn)只要計(jì)算一次。( )0001( , )( ,)()!kkkx yx

23、yyyk0( )00( , )( , ) ( , )1( , )( ,) ( , )()!ffkkkx y dyf x k g y kf x kx yg y kyydyk1.核函數(shù)的多極展開（平面問題常值元為例）*( )000001( )000010( )001( , )( )Re(, )(, )Re(, )(, ) Re(, )(, )Im(, )(, ) Im(, )(, )jjirirefirefrkkiirefiiimfkkriimfrkuTx y dyf xy k cy kfxy k cy kfxy k cy kf xy k cy kfxy k cy k( )000Im(, )(,

24、)iiimfikfxy k cy k( )( )ln( ),0Re( )Im( )( , ) ( , ), ( , ),1,2,.1,1,2,.rikkkx kxxf x kfx kfx kkxxkx基層的郵局可以方便地將任何地方轉(zhuǎn)來的郵件送給管轄區(qū)域任何住戶展成復(fù)數(shù)泰勒級數(shù)2.多極展開系數(shù)的轉(zhuǎn)移（M2M）*( )111101( )111110( )111( , )( )Re(, )(, )Re(, )(, ) Re(, )(, )Im(, )(, ) Im(, )(, )ImjjirirefirefrkkiirefiiimfkkrirefrkuTx y dyf xy k cy kfxy k

25、cy kfxy k cy kf xy k cy kfxy k cy k( )110(, )(, )iirefikfxy k cy k01100110100010001(, ),01(, )(,0)()(,)Re()(,)Im()(,) () ,1,2,.refkkqrefrefkrefrefrqqreficy kkcy kcyyyCcy kqyy cy kqkyy cy kqyyk 上級郵局利用下級郵局可以方便地將任何地方轉(zhuǎn)來的郵件送給管轄區(qū)域任何住戶3.多極展開系數(shù)向局部展開系數(shù)的轉(zhuǎn)移（M2L）*0000000000000000( , )( )Re ()(, )Re Re()()(, ) R

26、e Im()()(, )Im ()(, ) Re()()(, )Im Im()()(jjikkrenrenrkkkreniimnkkimnrimniuTx y dyxxdx kxxxxdx kxxxxdx kxxdx kxxxxdx kxxxxdx0, )k000000000001000101000000000001(,0)(,0)ln()(, )Re()(, )Im()(, )()(,0)11(, )(, )Re()(, )Im()()()prenrefrefrefrrefikkrefkrenl krefrefrllkdxcyxycy kxy cy kxy cy kxycydx lCcy k

27、xy cy kxl yxyxxy 001)(, ) 1,2,.prefiky cy kl發(fā)信方上級郵局利用收信方郵局可以將任何地方轉(zhuǎn)來的郵件送達(dá)收信的任何住戶4.局部展開系數(shù)的轉(zhuǎn)移（L2L）*1111101111111111( , )( )Re ()( , )Re Re()()( , ) Re Im()()( , )Im ()( , ) Re()()( , )Im Im()()(jjikkrenrenrkkkreniimnkkimnrimniuTx y dyxxdx kxxxxdx kxxxxdx kxxdx kxxxxdx kxxxxdx1, )k1100100100( , )()(, )R

28、e()(, )Im()(, ) 0,1,2,.plk lrenkrenrenrrenik ldx lCxxdx kxx dx kxx dx kl發(fā)信方將郵件交給最近的基層郵局就可以利用郵局系統(tǒng)送達(dá)任何收信人快速多極邊界元法采用迭代解法u 常規(guī)邊界元法解題規(guī)模較小，自由度數(shù)通常小于一萬，常用Gauss消去法解非對稱滿陣方程組。u 快速多極邊界元法，解題規(guī)模達(dá)到數(shù)百萬、上千萬，并不形成滿陣方程組，不存儲方程組系數(shù)，一定配合使用迭代解法。u 為了保證迭代盡快收斂，采用適當(dāng)?shù)念A(yù)條件處理技術(shù)也十分重要。快速多極邊界元法考題n計(jì)算考例（和邊界元法比較）多極快速邊界元法考題計(jì)算模型計(jì)算模型復(fù)數(shù)泰勒級數(shù)復(fù)數(shù)泰勒級數(shù)項(xiàng)數(shù)項(xiàng)數(shù)迭代次數(shù)迭代次數(shù)自由度數(shù)自由度數(shù)CPUCPU時間時間存儲需求存儲需求MB數(shù)數(shù)2 22 22525303026,88026,88000:19:5100:19:51505040404040252555555