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

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

2、效的工程數(shù)值分析方法在固體力學(xué)領(lǐng)域有限元法最重要的補(bǔ)充邊界元法間接法位勢問題（Smith & Pierce, 1958)彈性力學(xué)（Massonet，1965）邊界元法直接法位勢問題（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)用最多的部門也從未超過有限元法的十分之一,研究組邊界元研究歷史回顧,我們研究組邊界元法研究開始于1979年基于彈性力學(xué)問題Rizzo型邊界積分方程－邊界元法研究了彈性應(yīng)力集中問題和薄板彎曲問題研究了邊界元－有限元耦合方法研究了邊界元法在形狀優(yōu)化缺陷識別等逆問題中的應(yīng)用Galerkin對稱邊界元法用于結(jié)構(gòu)極限與安定分

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

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

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

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

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

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

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

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

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

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

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

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

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

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

18、,在復(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ī)分布不同形狀夾雜二維彈性體的變形圖,,

19、在復(fù)合材料模擬中的應(yīng)用,含100個隨機(jī)分布不同形狀夾雜二維彈性體的應(yīng)力分布,,在復(fù)合材料模擬中的應(yīng)用,含100個隨機(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)最初針對問題：多粒子相互作用的有勢場 存儲和計算復(fù)雜度：O(N2)Barnes, Hut, (1986) 樹結(jié)構(gòu)算法 存儲和計算復(fù)雜度：O(NlogN)Greengard, Rokhlin, (1986) 快速多極算法

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

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

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

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

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

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

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

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

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