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

1、邊界元法的若干進(jìn)展和它在固體力學(xué)中的應(yīng)用,清華大學(xué)工程力學(xué)系 姚振漢,引 言,彈性力學(xué)的三種提法,微分提法,變分提法,積分提法,偏微分方程邊值問(wèn)題,泛函極值問(wèn)題,邊界積分方程問(wèn)題,?,求解析解,差分法求數(shù)值解,李茲法求近似解,有限元法求數(shù)值解,邊界元法求數(shù)值解,?,?,歸結(jié)為求解線性代數(shù)方程組,?；癁槌Ｎ⒎址匠?,三種提法是完全等價(jià)的,邊界元發(fā)展歷史回顧,邊界積分方程－邊界元法有限元法（1955，56）之后發(fā)展起來(lái)的一種精確高

2、效的工程數(shù)值分析方法在固體力學(xué)領(lǐng)域有限元法最重要的補(bǔ)充邊界元法間接法位勢(shì)問(wèn)題（Smith & Pierce, 1958)彈性力學(xué)（Massonet，1965）邊界元法直接法位勢(shì)問(wèn)題（Jaswon，1963）彈性力學(xué)（Rizzo，1967）,邊界元發(fā)展歷史回顧,1994－2003被SCI收錄的論文與邊界元法有關(guān)的有3904篇與有限元法有關(guān)的為16823篇1990－2002被EI收錄的論文與邊界

3、元法有關(guān)的有19968篇與有限元法有關(guān)的為75184篇與斷裂力學(xué)有關(guān)的為23647篇在工程應(yīng)用方面在應(yīng)用最多的部門(mén)也從未超過(guò)有限元法的十分之一,研究組邊界元研究歷史回顧,我們研究組邊界元法研究開(kāi)始于1979年基于彈性力學(xué)問(wèn)題Rizzo型邊界積分方程－邊界元法研究了彈性應(yīng)力集中問(wèn)題和薄板彎曲問(wèn)題研究了邊界元－有限元耦合方法研究了邊界元法在形狀優(yōu)化缺陷識(shí)別等逆問(wèn)題中的應(yīng)用Galerkin對(duì)稱邊界元法用于結(jié)構(gòu)極限與安定分

4、析等問(wèn)題研究了精確高效的計(jì)算方法，提出了邊界元法誤差的一種直接估計(jì)。對(duì)于彈性接觸問(wèn)題，提出了單元與單元間協(xié)調(diào)的接觸方案，研究了二維、三維移動(dòng)、滾動(dòng)接觸。,研究組邊界元研究近年工作,2000前后針對(duì)復(fù)合材料，對(duì)于含隨機(jī)分布大量夾雜的二維彈性固體提出了一種重復(fù)相似子域邊界元法，計(jì)算了100多個(gè)夾雜、近萬(wàn)自由度的問(wèn)題。研究了在微機(jī)機(jī)群上的并行算法。2000年在由8臺(tái)微機(jī)組成的機(jī)群上最大計(jì)算規(guī)模45,000自由度。近年來(lái)多極快速算法在邊

5、界元法中的應(yīng)用給邊界元法解決復(fù)雜工程與科學(xué)問(wèn)題展示了廣闊的前景。用于含隨機(jī)分布夾雜二維、三維彈性體數(shù)值模擬，一臺(tái)微機(jī)可計(jì)算數(shù)十萬(wàn)自由度的問(wèn)題，在微機(jī)機(jī)群并行系統(tǒng)最大的二維算例有800萬(wàn)自由度。計(jì)算了含16384條隨機(jī)分布裂紋的二維無(wú)限彈性體，1,572,864自由度，研究了相應(yīng)的裂紋擴(kuò)展問(wèn)題。,鑒于邊界元法始終是計(jì)算力學(xué)中一個(gè)重要的研究領(lǐng)域快速多極邊界元法是近年受到特別關(guān)注的一個(gè)研究方向本人已經(jīng)從事邊界元法研究25年，幾乎沒(méi)有

6、間斷近年研究組研究重點(diǎn)放在快速多極邊界元法的研究方面因此今天就以此為題，希望能對(duì)各位同學(xué)的學(xué)習(xí)有一定的幫助。,彈性力學(xué)邊界積分方程,彈性靜力學(xué)的邊界積分方程其中,,,,,,,,,,,對(duì)于同一彈性體的兩種變形狀態(tài),,?,?,,,,通常一個(gè)狀態(tài)是待求狀態(tài)另一狀態(tài)是已知的輔助狀態(tài),彈性力學(xué)的邊界積分方程可以由Betti定理出發(fā)導(dǎo)出,,,,,以無(wú)限彈性空間中一點(diǎn)作用單位集中力的Kelvin問(wèn)題的解為輔助解,,?,,由Bett

7、i定理可得Somigliana等式（無(wú)體積力情況）,,,再將集中力作用點(diǎn)（源點(diǎn)）P 趨于邊界點(diǎn) p，考慮到源點(diǎn)是積分的奇異點(diǎn)，作適當(dāng)處理，即可得到前面給出的邊界積分方程。,彈性力學(xué)邊界元法,將邊界分成邊界單元在每個(gè)單元上將邊界變量插值離散可得,,,,將弱解代入方程得到誤差對(duì)于加權(quán)余量法的配點(diǎn)格式，權(quán)函數(shù)采用Dirac-delta函數(shù)，要求即得,,,,彈性力學(xué)邊界元法,將方程寫(xiě)成矩陣形式，得其中H，G矩

8、陣元素由核函數(shù)與形函數(shù)在單元上的積分求得。將邊界變量列矢量 U，T 按未知量與給定量重新排列，可得邊界元法的求解代數(shù)方程組即,,,,,,,彈性力學(xué)邊界元法,核函數(shù)與形函數(shù)乘積在單元上的積分矩陣元素都是核函數(shù)與形函數(shù)乘積在單元上的積分，矩陣是滿陣。主要計(jì)算量就是計(jì)算這些積分，以及求解滿陣代數(shù)方程組。對(duì)于規(guī)模不太大的問(wèn)題，計(jì)算積分的工作量是主要計(jì)算量。非奇異積分采用等精度高斯積分格式求積，高斯點(diǎn)數(shù)由精度要求對(duì)不同情

9、況自動(dòng)確定。當(dāng)核函數(shù)的源點(diǎn)落入積分單元時(shí)出現(xiàn)奇異積分，包括弱奇異積分和柯西主值積分。,彈性力學(xué)邊界元法,邊界元法的優(yōu)缺點(diǎn),邊界元法與有限元法及其它數(shù)值方法相比較的優(yōu)缺點(diǎn)優(yōu)點(diǎn)：高精度（由于采用了解析基本解）降低了維數(shù)，便于模擬復(fù)雜邊界形狀對(duì)于高梯度、甚至有奇異性的問(wèn)題，不僅有較高的精度，而且在同等精度條件下有較高效率適合于處理無(wú)限域、半無(wú)限域問(wèn)題適合于處理彈性接觸等邊界條件非線性問(wèn)題缺點(diǎn)：適用范圍遠(yuǎn)沒(méi)有有限元法廣泛

10、解題規(guī)模受限制（方程系數(shù)矩陣為滿陣）對(duì)于域內(nèi)方程非線性問(wèn)題優(yōu)勢(shì)減弱或喪失,邊界元法的特點(diǎn),早期有的文章通過(guò)簡(jiǎn)單算例高估了邊界元法的計(jì)算精度，聲稱能達(dá)到0.01%，甚至更高。這里要分兩類問(wèn)題，一類是沒(méi)有離散誤差的問(wèn)題，稱為簡(jiǎn)單問(wèn)題，離散插值的邊界變量能精確滿足給定邊界條件，此時(shí)能達(dá)到很高精度，例如10－6，甚至10－8。對(duì)于一般問(wèn)題，有離散誤差，只有合理劃分足夠多的邊界單元，才能達(dá)到要求的計(jì)算精度。薄板梁純彎曲問(wèn)題，只要

11、采用二次單元，就是一個(gè)簡(jiǎn)單問(wèn)題。通常認(rèn)為細(xì)長(zhǎng)薄板梁要?jiǎng)澐州^多單元才能達(dá)到滿意精度，我們只用4個(gè)二次元，就達(dá)10－4精度。,有限元法中單元邊長(zhǎng)比不能太大，相鄰單元尺寸也不能相差太大。邊界元法則并不受此限制，上述100:1純彎薄板梁相鄰單元長(zhǎng)度比為100，只要保證積分精度等運(yùn)算精度，還是可以得到高精度的結(jié)果。有限元法中高斯積分通常只用1－2，或2?2個(gè)高斯積分點(diǎn)。邊界元法常采用等精度高斯積分，根據(jù)給定積分精度要求來(lái)確定用多少高

12、斯點(diǎn)。對(duì)于上述細(xì)長(zhǎng)薄板梁的算例，源點(diǎn)在邊長(zhǎng)為1的短邊中點(diǎn)、在長(zhǎng)邊二次單元上的積分，需要50－60個(gè)高斯點(diǎn)進(jìn)行高斯積分。,邊界元法的特點(diǎn),邊界元法在解出未知的邊界變量之后，通過(guò)域內(nèi)變量的邊界積分公式可以求得域內(nèi)任意點(diǎn)的位移、應(yīng)變和應(yīng)力。這樣得到的域內(nèi)位移場(chǎng)精確滿足Navier方程，但是當(dāng)它趨于邊界時(shí)的極限一般說(shuō)來(lái)和邊界給定量或邊界元法解出的邊界變量并不一致。這種差別是和離散誤差密切相關(guān)的。對(duì)于沒(méi)有離散誤差的簡(jiǎn)單問(wèn)題，域內(nèi)解趨于邊界的極

13、限和邊界變量是一致的，只有很小的運(yùn)算誤差。這也導(dǎo)致在有離散誤差的情況下，在近邊界區(qū)求得的內(nèi)點(diǎn)變量有較大的誤差，即通常稱為“邊界效應(yīng)”。,邊界元法的特點(diǎn),邊界元法模擬光彈實(shí)驗(yàn),為了便于將得到的結(jié)果和工程界熟悉的光彈性實(shí)驗(yàn)結(jié)果比較，域內(nèi)應(yīng)力分布也畫(huà)出了等色線和等傾線圖。而且我們一再表明，在數(shù)值模擬已經(jīng)很成熟的今天，進(jìn)行光彈性實(shí)驗(yàn)這樣的物理模擬已經(jīng)沒(méi)有多大意義。因?yàn)樗膊皇菍?shí)測(cè)，光彈性模型和實(shí)際問(wèn)題相比引進(jìn)了許多誤差，得到的結(jié)果并不比

14、數(shù)值模擬結(jié)果更接近真實(shí)情況；而且數(shù)值模擬比光彈性模擬更經(jīng)濟(jì)便捷。下圖是帶孔板單向拉伸的等色線圖（條紋等級(jí)反映主應(yīng)力差）。,邊界元法模擬光彈實(shí)驗(yàn),邊界元法并行計(jì)算的工程應(yīng)用,邊界元法結(jié)合并行計(jì)算，在8臺(tái)微機(jī)組成的機(jī)群并行系統(tǒng)上計(jì)算了一些規(guī)模稍大的工程實(shí)際問(wèn)題。例如：石油鉆桿的偏心鉆挺，49818自由度，8臺(tái)微機(jī)，計(jì)算43小時(shí)41分鐘。渤海石油公司“濱海109”海洋鋪管船的錨機(jī)滾筒軸，24444自由度，7臺(tái)微機(jī)，計(jì)算436分鐘。

15、這是2000年博士生尹欣得到的結(jié)果，當(dāng)時(shí)文獻(xiàn)中最大計(jì)算規(guī)模10萬(wàn)自由度。,子域 1的網(wǎng)格劃分,子域 2~5的網(wǎng)格劃分,偏心鉆挺模型（5萬(wàn)自由度）,有限元模型（16000結(jié)點(diǎn)，13000千單元）,鉆挺頂部Mises應(yīng)力圖（邊界元法，單位：MPa）,鉆挺頂部Mises應(yīng)力圖（邊界元法，單位：MPa）,鉆挺頂部Mises應(yīng)力圖（有限元法，單位：MPa）,鉆挺頂部Mises應(yīng)力圖（有限元法，單位：MPa）,滾筒軸邊界元模型：5個(gè)子域，8184結(jié)

16、點(diǎn)，1721單元,滾筒軸有限元模型：13000結(jié)點(diǎn)，64000單元,邊界元解：左端局部Mises應(yīng)力（單位：MPa）,有限元解：左端局部Mises應(yīng)力（單位：MPa）,邊界元解：右端局部Mises應(yīng)力（單位：MPa）,有限元解：右端局部Mises應(yīng)力（單位：MPa）,在復(fù)合材料模擬中的應(yīng)用,針對(duì)復(fù)合材料數(shù)值模擬，研究組對(duì)于含隨機(jī)分布孔洞或夾雜的二維彈性體進(jìn)行了深入的研究。提出了重復(fù)相似子域邊界元法，提高了計(jì)算效率。分析計(jì)算了含100個(gè)圓

17、孔、圓形夾雜、橢圓形夾雜、多種形狀?yuàn)A雜的二維彈性體，得到了等效材料常數(shù)和體積比、夾雜與基底材料常數(shù)比等的關(guān)系?？紤]了夾雜和基底之間有界面層的情況。對(duì)此類問(wèn)題，引進(jìn)快速算法之后，使計(jì)算規(guī)模大大提高，初步工作已經(jīng)從100多個(gè)夾雜增加到了1600個(gè)夾雜。,,在復(fù)合材料模擬中的應(yīng)用,重復(fù)相似子域邊界元法,,在復(fù)合材料模擬中的應(yīng)用,含100個(gè)隨機(jī)分布圓形夾雜的二維彈性體,,在復(fù)合材料模擬中的應(yīng)用,含100個(gè)隨機(jī)分布圓形夾雜的二維彈性體變形圖,

18、,在復(fù)合材料模擬中的應(yīng)用,含100個(gè)隨機(jī)分布圓形夾雜的二維彈性體應(yīng)力分布,,在復(fù)合材料模擬中的應(yīng)用,含100個(gè)隨機(jī)分布橢圓形夾雜的二維彈性體,,在復(fù)合材料模擬中的應(yīng)用,含100個(gè)隨機(jī)分布橢圓形夾雜二維彈性體的變形圖,,在復(fù)合材料模擬中的應(yīng)用,含100個(gè)隨機(jī)分布橢圓形夾雜二維彈性體的應(yīng)力分布,,在復(fù)合材料模擬中的應(yīng)用,含100個(gè)隨機(jī)分布不同形狀?yuàn)A雜的二維彈性體,,在復(fù)合材料模擬中的應(yīng)用,含100個(gè)隨機(jī)分布不同形狀?yuàn)A雜二維彈性體的變形圖,,

19、在復(fù)合材料模擬中的應(yīng)用,含100個(gè)隨機(jī)分布不同形狀?yuàn)A雜二維彈性體的應(yīng)力分布,,在復(fù)合材料模擬中的應(yīng)用,含100個(gè)隨機(jī)分布圓形夾雜（有界面層）的二維彈性體,,在復(fù)合材料模擬中的應(yīng)用,含圓形夾雜二維彈性體等效彈性模量和模量比的關(guān)系,,在復(fù)合材料模擬中的應(yīng)用,含圓形夾雜二維彈性體等效體積模量和模量比的關(guān)系,,在復(fù)合材料模擬中的應(yīng)用,含圓形夾雜二維彈性體等效剪切模量和模量比的關(guān)系,,在復(fù)合材料模擬中的應(yīng)用,含圓形夾雜二維彈性體等效彈性模量和模量

20、比的關(guān)系,快速多極算法,快速多極算法 (fast multipole method, FMM)最初針對(duì)問(wèn)題：多粒子相互作用的有勢(shì)場(chǎng) 存儲(chǔ)和計(jì)算復(fù)雜度：O(N2)Barnes, Hut, (1986) 樹(shù)結(jié)構(gòu)算法 存儲(chǔ)和計(jì)算復(fù)雜度：O(NlogN)Greengard, Rokhlin, (1986) 快速多極算法

21、 存儲(chǔ)和計(jì)算復(fù)雜度：O(N)Hrycak T, Rokhlin V, (1998) 新型快速多極算法 存儲(chǔ)和計(jì)算復(fù)雜度：O(N),快速多極邊界元法,常規(guī)邊界元法由于方程組的系數(shù)矩陣為非對(duì)稱滿陣，對(duì)于N自由度系統(tǒng)，其運(yùn)算量為N 3量級(jí)，對(duì)于存儲(chǔ)的要求為N 2量級(jí)。邊界元快速算法通過(guò)快速多極展開(kāi)，使運(yùn)算量和存儲(chǔ)量均減少到NlogN 量級(jí)。最近引入了局部展開(kāi)的思想，進(jìn)一步使運(yùn)算量和存儲(chǔ)

22、量要求減少到了N 量級(jí)。研究組初步的工作已經(jīng)把解題規(guī)模明顯擴(kuò)大，由原來(lái)8臺(tái)微機(jī)并行最多計(jì)算49,818自由度問(wèn)題，到用邊界元快速算法一臺(tái)微機(jī)計(jì)算544,000自由度和1,572,864自由度。,邊界元快速算法和常規(guī)邊界元法的比較,快速多極邊界元法基本思想,快速多極算法的基本思想：如果每個(gè)人要自己把信送到每個(gè)收信人，是非常費(fèi)時(shí)的，當(dāng)收信人的范圍很大時(shí)是不可能的。但是現(xiàn)代的郵政系統(tǒng)，使每個(gè)人可以方便地給世界各地的收信人發(fā)信。

23、傳統(tǒng)邊界元法形成滿陣方程組的系數(shù)要一個(gè)個(gè)地獨(dú)立計(jì)算，就像每個(gè)人親自送信一樣?？焖俣鄻O算法就像現(xiàn)代郵政系統(tǒng)一樣大大提高了效率。,,樹(shù)型存儲(chǔ)結(jié)構(gòu),將求解區(qū)域與邊界有關(guān)部分劃分為四叉樹(shù)（二維問(wèn)題）或八叉樹(shù)（三維問(wèn)題）,,,近場(chǎng)與遠(yuǎn)場(chǎng),對(duì)源點(diǎn) x，近場(chǎng)按常規(guī)邊界元法計(jì)算方程系數(shù)，遠(yuǎn)場(chǎng)用快速多極展開(kāi)，不再計(jì)算一個(gè)個(gè)系數(shù)。,,,快速多極算法基本步驟,對(duì)于遠(yuǎn)場(chǎng)的快速多極算法有四個(gè)基本步驟：,,核函數(shù)的多極展開(kāi) 多極展開(kāi)系數(shù)（多極

24、矩）g (y, k) 對(duì)所有核函數(shù)源點(diǎn)只要計(jì)算一次。,,,,核函數(shù)的多極展開(kāi)（平面問(wèn)題常值元為例）,,,,,,,基層的郵局可以方便地將任何地方轉(zhuǎn)來(lái)的郵件送給管轄區(qū)域任何住戶,展成復(fù)數(shù)泰勒級(jí)數(shù),,多極展開(kāi)系數(shù)的轉(zhuǎn)移（M2M）,,,,,,,,上級(jí)郵局利用下級(jí)郵局可以方便地將任何地方轉(zhuǎn)來(lái)的郵件送給管轄區(qū)域任何住戶,多極展開(kāi)系數(shù)向局部展開(kāi)系數(shù)的轉(zhuǎn)移（M2L）,,,,,,,,,,,發(fā)信方上級(jí)郵局利用收信方郵局可以將任何地方轉(zhuǎn)來(lái)的郵件送達(dá)收信的任何

25、住戶,局部展開(kāi)系數(shù)的轉(zhuǎn)移（L2L）,,,,,,,,,,,發(fā)信方將郵件交給最近的基層郵局就可以利用郵局系統(tǒng)送達(dá)任何收信人,快速多極邊界元法采用迭代解法,常規(guī)邊界元法解題規(guī)模較小，自由度數(shù)通常小于一萬(wàn)，常用Gauss消去法解非對(duì)稱滿陣方程組。,快速多極邊界元法，解題規(guī)模達(dá)到數(shù)百萬(wàn)、上千萬(wàn)，并不形成滿陣方程組，不存儲(chǔ)方程組系數(shù)，一定配合使用迭代解法。,為了保證迭代盡快收斂，采用適當(dāng)?shù)念A(yù)條件處理技術(shù)也十分重要。,快速多極邊界元法考題,計(jì)算考例（

26、和邊界元法比較）,,多極快速邊界元法考題,,,快速多極邊界元法算例,含1600個(gè)隨機(jī)分別圓形夾雜的二維彈性體,,Eb =200.0MPa,快速多極邊界元法算例,夾雜體積比、模量比對(duì)等效體積模量的影響,,快速多極邊界元法算例,快速算法老方案和常規(guī)邊界元法計(jì)算時(shí)間比較,,快速多極邊界元法算例,快速算法老方案和常規(guī)邊界元法存儲(chǔ)需求比較,,快速多極邊界元法算例,快速算法新方案和老方案計(jì)算時(shí)間比較,,快速多極邊界元法算例,含300個(gè)隨機(jī)分布球形夾

27、雜的三維彈性體,,快速多極邊界元法算例,300個(gè)隨機(jī)分布球形夾雜的界面正應(yīng)力,,快速多極邊界元法算例,300個(gè)隨機(jī)分布球形夾雜的界面剪應(yīng)力,,快速多極邊界元法算例,等效彈性模量與夾雜模量比的關(guān)系,,快速多極邊界元法算例,等效體積模量與夾雜模量比的關(guān)系,,快速多極邊界元法算例,等效彈性模量與夾雜模量比的關(guān)系,,快速多極邊界元法算例,等效體積模量與夾雜模量比的關(guān)系,,快速多極邊界元法算例,含100根隨機(jī)分布同方向長(zhǎng)纖維三維彈性體（長(zhǎng)細(xì)比＝1

28、1）,,快速多極邊界元法算例,100根隨機(jī)分布同方向長(zhǎng)纖維界面正應(yīng)力分布,,快速多極邊界元法算例,100根隨機(jī)分布同方向長(zhǎng)纖維界面剪應(yīng)力分布,,快速多極邊界元法算例,含100根隨機(jī)分布長(zhǎng)纖維三維彈性體（最大方位偏差30?和90?）,,快速多極邊界元法算例,含300根隨機(jī)分布短纖維三維彈性體（方位偏差30?）,,快速多極邊界元法算例,300根隨機(jī)分布短纖維界面剪應(yīng)力分布,,快速多極邊界元法算例,等效彈性模量與纖維長(zhǎng)細(xì)比、方位角偏差的關(guān)系,

29、,快速多極邊界元法算例,等效體積模量與纖維長(zhǎng)細(xì)比、方位角偏差的關(guān)系,,快速多極邊界元法算例,含10000條隨機(jī)分布裂紋的二維彈性體雙向均勻拉伸,,快速多極邊界元法算例,左下角區(qū)域裂紋張開(kāi)位移,快速多極邊界元法算例,,含100條裂紋方板中裂紋的擴(kuò)展,,快速多極邊界元法算例,含100條裂紋方板中裂紋的擴(kuò)展,,快速多極邊界元法算例,含100條裂紋方板中裂紋的擴(kuò)展,,快速多極邊界元法算例,含100條裂紋方板中裂紋的擴(kuò)展,研究展望,多極快速邊界元