第1章緒論
1.1數(shù)值方法概述隨著科學(xué)技術(shù)的不斷發(fā)展，需要解決的工程問(wèn)題也越來(lái)越復(fù)雜，對(duì)于大多數(shù)問(wèn)題，由于求解問(wèn)題幾何形狀的復(fù)雜性或計(jì)算介質(zhì)的非線性，人們已經(jīng)很難得到問(wèn)題的解析答案。另外，隨著計(jì)算機(jī)性能的日益提高，求解工程問(wèn)題的數(shù)值計(jì)算方法也不斷成熟，現(xiàn)在幾乎所有的大型工程問(wèn)題都要借助數(shù)值計(jì)算進(jìn)行分析或評(píng)估，為科技人員的工程設(shè)計(jì)提供依據(jù)或參考。
目前，已經(jīng)發(fā)展起來(lái)的用得較多的數(shù)值計(jì)算方法有五大類(lèi)：有限差分法、有限體積法、有限單元法、無(wú)網(wǎng)格法和邊界元法。
(1) 有限差分法(finite difference method)［1］。有限差分法是將所考慮的區(qū)域劃分成網(wǎng)格，用差分近似微分，把微分方程變成差分方程。也就是通過(guò)數(shù)學(xué)上的近似，把求解微分方程的問(wèn)題變換成求解關(guān)于節(jié)點(diǎn)未知量的代數(shù)方程問(wèn)題。該方法簡(jiǎn)單、易懂，便于復(fù)雜微分方程的求解，因此在流體力學(xué)領(lǐng)域被廣泛采用［2］。但當(dāng)求解問(wèn)題的幾何形狀復(fù)雜時(shí)，按空間坐標(biāo)的離散變得困難，網(wǎng)格的正交性不易保留，導(dǎo)致計(jì)算精度降低。
(2) 有限體積法(finite volume method)［3］。有限體積法是基于物理問(wèn)題控制方程的積分形式的數(shù)值方法。解題思路是：把計(jì)算域離散成有限個(gè)互不重疊的控制體單元（網(wǎng)格），通過(guò)將積分形式的控制方程作用于每一個(gè)體單元來(lái)建立離散的代數(shù)方程組。有限體積法既有有限差分法的特點(diǎn)，又有有限單元法的特點(diǎn)，方法簡(jiǎn)單、幾何適應(yīng)性好，是現(xiàn)代計(jì)算流體力學(xué)中占主導(dǎo)地位的數(shù)值方法［4］。其缺點(diǎn)是物理量的空間導(dǎo)數(shù)相關(guān)量（如通量）是由周?chē)�體單元中心處的值決定的，因此精度較低，特別是靠近邊界的通量值更難計(jì)算準(zhǔn)確。
(3) 有限單元法(finite element method)［5］。有限單元法是通過(guò)變分原理建立含權(quán)函數(shù)（形函數(shù)）的體積分形式方程的方法。解題思路是：首先將計(jì)算域離散成有限個(gè)互不重疊的有一定規(guī)則的節(jié)點(diǎn)組成的單元，然后對(duì)每個(gè)單元積分并通過(guò)組集形成總體代數(shù)方程組。有限單元法的單元可以按不同的連接方式進(jìn)行組合，每個(gè)單元可以有不同的形狀和材料性質(zhì)，因此有限單元法具有幾何適應(yīng)性強(qiáng)、可靈活處理不同物性參數(shù)的優(yōu)點(diǎn)，在各個(gè)領(lǐng)域都得到了廣泛應(yīng)用［6］。有限單元法的缺點(diǎn)是：物理量的空間導(dǎo)數(shù)是通過(guò)對(duì)形函數(shù)的求導(dǎo)得到的，因此精度要比物理量本身低一階；在金屬成形、優(yōu)化計(jì)算、滲流問(wèn)題自由面的確定等運(yùn)動(dòng)邊界問(wèn)題中，有限元網(wǎng)格可能產(chǎn)生畸變和重疊，以致計(jì)算精度下降或計(jì)算中止。
(4) 無(wú)網(wǎng)格法(meshless method)［7］。無(wú)網(wǎng)格法是基于構(gòu)造點(diǎn)插值函數(shù)的數(shù)值方法，因此只需要在計(jì)算域內(nèi)布置一系列的離散點(diǎn)即可，不需要網(wǎng)格單元，具有很強(qiáng)的解題靈活性［8，9］。但無(wú)網(wǎng)格法發(fā)展得還不夠成熟，缺少堅(jiān)實(shí)的理論基礎(chǔ)和嚴(yán)格的數(shù)學(xué)證明，因此計(jì)算精度、守恒性等一直沒(méi)有明確的答案。此外，無(wú)網(wǎng)格法是基于點(diǎn)的算法，因此布點(diǎn)數(shù)量和方案會(huì)影響計(jì)算時(shí)間和精度。另外，由于不使用單元，對(duì)于幾何較復(fù)雜的運(yùn)動(dòng)邊界問(wèn)題，邊界變化時(shí)要判斷重新分布后的點(diǎn)是內(nèi)部點(diǎn)、外部點(diǎn)還是邊界點(diǎn)，這時(shí)會(huì)存在困難。
(5) 邊界元法(boundary element method，BEM)［10］。邊界元法是基于格林公式和問(wèn)題的基本解將控制微分方程轉(zhuǎn)化為邊界積分方程的一種數(shù)值方法。其主要優(yōu)點(diǎn)是：①只需要將邊界離散成單元，因而準(zhǔn)備數(shù)據(jù)簡(jiǎn)單、便于復(fù)雜幾何問(wèn)題的建模［11］；②能夠自動(dòng)滿(mǎn)足無(wú)限遠(yuǎn)處的邊界條件，因而適合于求解無(wú)限與半無(wú)限域問(wèn)題［12］；③所求物理量的空間導(dǎo)數(shù)的計(jì)算公式可以解析地從基本邊界積分方程中導(dǎo)出， 因此所求與導(dǎo)數(shù)相關(guān)的物理量（如通量、應(yīng)力等）與物理量本身具有同樣級(jí)別的精度［13］；④在求解運(yùn)動(dòng)邊界問(wèn)題時(shí)，移動(dòng)邊界節(jié)點(diǎn)的位移與原邊界節(jié)點(diǎn)的坐標(biāo)相加就自然形成了新的邊界單元信息，不需要專(zhuān)門(mén)重構(gòu)單元，也不會(huì)有網(wǎng)格重疊的問(wèn)題［14��16］。此外，由于在計(jì)算域的邊界上有單元信息可用，所以通過(guò)單元積分很容易判斷一點(diǎn)是內(nèi)部、外部還是邊界點(diǎn)。基于這些優(yōu)點(diǎn)，邊界元法在一些領(lǐng)域（如運(yùn)動(dòng)邊界［15��17］、裂紋［18，19］、接觸［20，21］、輻射［22，23］、超薄結(jié)構(gòu)［24��28］、無(wú)限域和半無(wú)限域［14，29，30］、聲學(xué)［31］等）的應(yīng)用優(yōu)于有限元法。
1.2邊界元法的發(fā)展歷史
邊界元法已發(fā)展成為求解工程與科學(xué)問(wèn)題的常用數(shù)值分析方法之一。它是一種在經(jīng)典的積分方程基礎(chǔ)上，吸收了有限元法的離散技術(shù)而發(fā)展起來(lái)的數(shù)值方法。邊界元法通過(guò)采用一個(gè)滿(mǎn)足無(wú)限或半無(wú)限域場(chǎng)方程的奇異函數(shù)——基本解作為權(quán)函數(shù)，將基于問(wèn)題控制方程的域積分方程轉(zhuǎn)化為邊界積分方程，并將邊界離散成邊界單元來(lái)求解邊界未知量的數(shù)值解。
邊界元法的產(chǎn)生可追溯到19世紀(jì)，當(dāng)時(shí)有人提出了積分等式和位勢(shì)理論的數(shù)學(xué)問(wèn)題，把線性偏微分方程的邊值問(wèn)題轉(zhuǎn)化為邊界積分方程求解。1905年，F(xiàn)redholm將積分方程應(yīng)用于求解彈性力學(xué)問(wèn)題［32］。1953年，Muskhelishvili［33］和Kellogg［34］分別將積分方程法用于求解結(jié)構(gòu)力學(xué)和位勢(shì)問(wèn)題。1965年，Mikhlin［35］解決了積分方程理論中的奇異性問(wèn)題，為積分方程法在工程中的應(yīng)用開(kāi)辟了道路。
20世紀(jì)60年代，電子計(jì)算機(jī)的發(fā)展開(kāi)創(chuàng)了數(shù)值求解積分方程的新時(shí)代，積分方程法作為數(shù)值計(jì)算方法開(kāi)始應(yīng)用于實(shí)際問(wèn)題。1963年，Jaswon［36］采用間接邊界積分方程法，成功地求解了位勢(shì)問(wèn)題和彈性力學(xué)問(wèn)題。這種方法的主要思想是沿邊界配置一種虛設(shè)的點(diǎn)源密度函數(shù)，先確定密度函數(shù)，再求邊界上的未知物理量。其缺點(diǎn)是待求的點(diǎn)源分布函數(shù)是虛構(gòu)的，不具有明確的物理意義，因此求解過(guò)程需要兩步完成。但它具有場(chǎng)量方程和場(chǎng)量梯度方程相互獨(dú)立的優(yōu)點(diǎn)，因而易于組合求解各種復(fù)雜邊界條件的邊值問(wèn)題。60年代，一些蘇聯(lián)學(xué)者對(duì)積分方程尤其是奇異積分方程的理論作了更為深入的研究，為進(jìn)一步應(yīng)用邊界積分方程方法開(kāi)辟了道路。與此同時(shí)，高速大型計(jì)算機(jī)的出現(xiàn)及其硬件的迅猛發(fā)展使離散求解積分方程成為可能。到了60年代后半期，邊界元法的研究更趨活躍，邊界元法的直接法被應(yīng)用于求解各類(lèi)問(wèn)題。1967年，Rizzo［37］用直接邊界元法求解了二維彈性問(wèn)題。1969年，Cruse［38］將此法推廣到三維彈性力學(xué)問(wèn)題。在直接法中，表述邊界積分方程的未知量是真實(shí)的物理量，通過(guò)求解積分方程可以直接得到邊界上所求的未知物理量。1975年，Cruse和Rizzo［39］出版了第一部邊界積分方程法的著作。1977年，Banerjee和Butterfield［40］首次采用了boundary element method這一名稱(chēng)。1978年，Brebbia在英國(guó)南安普頓召開(kāi)了第一屆國(guó)際邊界元法會(huì)議，出版了專(zhuān)著The Boundary Element Method for Engineers［41］。書(shū)中用加權(quán)余量法推導(dǎo)出了邊界積分方程，并指出加權(quán)余量法是最普遍的數(shù)值方法，如果以開(kāi)爾文（Kelvin）解作為權(quán)函數(shù)，從加權(quán)余量法可導(dǎo)出彈性力學(xué)問(wèn)題的邊界積分方程，通過(guò)將邊界離散成單元的方法可數(shù)值求解積分方程。至此，邊界元法這一名稱(chēng)得到了國(guó)際公認(rèn)。
自1978年第一屆國(guó)際邊界元法會(huì)議后，邊界元法會(huì)議幾乎每年在世界各地舉辦。世界各國(guó)已從基本理論與方法的研究向深廣領(lǐng)域發(fā)展，大量論文和專(zhuān)著先后問(wèn)世。在此時(shí)期，邊界元法在數(shù)學(xué)分析理論和數(shù)值積分方法的研究、邊界元法的完善和應(yīng)用范圍的拓寬以及邊界元法應(yīng)用軟件等方面均得到飛速發(fā)展。邊界元法的應(yīng)用遍及固體力學(xué)、流體力學(xué)、波動(dòng)學(xué)、傳熱學(xué)、電磁學(xué)等學(xué)科領(lǐng)域。
我國(guó)關(guān)于邊界元法的研究大約開(kāi)始于1978年，當(dāng)時(shí)杜慶華在國(guó)內(nèi)首先開(kāi)創(chuàng)了工程中邊界元法的研究，開(kāi)始跟蹤國(guó)際上這一領(lǐng)域的最新進(jìn)展，其研究領(lǐng)域主要在固體力學(xué)方面［42］。與此同時(shí)，王泳嘉［43］、鄭穎人等［44］開(kāi)始了巖土工程中的邊界元法研究，楊德全和趙忠生［45］在流體力學(xué)邊界元法方面開(kāi)了先河。值得一提的是，岑章志［46］對(duì)我國(guó)的彈塑性邊界元法做了開(kāi)拓性的研究工作，高效偉和Davies［13］發(fā)表了國(guó)際上第一個(gè)彈塑性邊界元程序，結(jié)束了三十多年來(lái)在非線性力學(xué)邊界元分析方面只有論文發(fā)表、沒(méi)有程序公布的局面。
近年來(lái)，我國(guó)學(xué)者在快速多極邊界元法研究方面取得了一系列的研究成果，代表性工作可見(jiàn)姚振漢與王海濤的著作［11］。此外，我國(guó)學(xué)者在近奇異積分計(jì)算［24��28］、多種介質(zhì)問(wèn)題［47��49］以及與CAD技術(shù)結(jié)合解決工程問(wèn)題［50］等方面的研究工作也引人注目。邊界元法的研究目前在我國(guó)正處于上升階段［51］，近年來(lái)在一些回國(guó)學(xué)者的帶領(lǐng)下越來(lái)越多的專(zhuān)家學(xué)者投身于邊界元法的研究中。
1.3邊界元法中的難點(diǎn)問(wèn)題及其研究進(jìn)展
邊界元法在解決接觸、斷裂力學(xué)、運(yùn)動(dòng)邊界、無(wú)限域與半無(wú)限域以及超薄結(jié)構(gòu)等問(wèn)題中具有獨(dú)特的優(yōu)勢(shì)，被大量地用于科學(xué)與工程問(wèn)題的計(jì)算分析，在許多應(yīng)用領(lǐng)域，邊界元法被公認(rèn)為有限元法的一個(gè)重要補(bǔ)充。然而，傳統(tǒng)的邊界元法有下述幾方面的弱點(diǎn)：
（1） 奇異性問(wèn)題。邊界元法中所用的基本解是奇異函數(shù)，數(shù)值計(jì)算時(shí)首先需要消除積分奇異性才能得到精確的計(jì)算結(jié)果。多年來(lái)已有大量的文獻(xiàn)在計(jì)算效率與精度方面提出了不少計(jì)算奇異積分的方法，如線性單元的解析消除法［10，52］、弱奇異積分的單元子分法［13，53，54］、強(qiáng)奇異積分的間接計(jì)算法［13，55］、高階奇異積分的直接計(jì)算法［56��58］等。這些方法在計(jì)算弱奇異和強(qiáng)奇異積分時(shí)非常有效，但在處理高階奇異積分時(shí)的穩(wěn)定性還需要進(jìn)一步考查。
（2） 滿(mǎn)系數(shù)矩陣問(wèn)題。邊界元法在求解問(wèn)題時(shí)所形成的方程組的系數(shù)矩陣是滿(mǎn)陣，因而占有較大的計(jì)算機(jī)內(nèi)存，難以解決大型工程問(wèn)題。為了解決此問(wèn)題，研究人員提出了兩種有效的解決方法，一種是快速多極法，另一種是區(qū)域分解法�？焖俣鄻O邊界元法［11，59］，通過(guò)將奇異核函數(shù)進(jìn)行級(jí)數(shù)展開(kāi)的技術(shù)，降低矩陣向量相乘操作的計(jì)算量級(jí)和存儲(chǔ)量級(jí)，達(dá)到節(jié)省內(nèi)存和提高計(jì)算速度的目的。區(qū)域分解法［10，60］的基本思想是把總求解域劃分成多個(gè)子域，對(duì)每個(gè)子域建立邊界元矩陣方程，然后利用子域間公共節(jié)點(diǎn)上的面力平衡條件和位移相容性條件把子域方程組集成總體系統(tǒng)方程組。這樣形成的系數(shù)矩陣是塊狀稀疏矩陣，可利用現(xiàn)有成熟的稀疏矩陣求解器（如LU分解法和GMRES迭代法（廣義最小殘量法））對(duì)系統(tǒng)方程組進(jìn)行有效求解。區(qū)域分解法是邊界元法中被廣泛應(yīng)用的方法，不僅能解決滿(mǎn)系數(shù)矩陣問(wèn)題，而且能通過(guò)按照材料性質(zhì)劃分子域的手段求解由不同材料組成的復(fù)合介質(zhì)問(wèn)題［61］，還能通過(guò)沿裂紋面劃分子域的方法求解斷裂力學(xué)問(wèn)題［19］。最近，作者基于區(qū)域分解法，提出了求解非均勻介質(zhì)問(wèn)題的三步變量凝聚技術(shù)［49］，形成的系統(tǒng)方程組只有公共節(jié)點(diǎn)位移作為未知量，可以求解大型工程問(wèn)題。雖然這些方程組的組建與求解技術(shù)能解決相當(dāng)規(guī)模的工程問(wèn)題，但對(duì)于超大型問(wèn)題，系統(tǒng)方程組的求解仍然是一個(gè)極具挑戰(zhàn)性的問(wèn)題。最近，高效偉等提出了求解非對(duì)稱(chēng)稀疏方程組的同時(shí)消元回代法求解技術(shù)［62，63］，在計(jì)算效率和儲(chǔ)存空間方面都有了顯著的提高。
（3） 非線性和非均質(zhì)問(wèn)題中的域積分問(wèn)題。傳統(tǒng)的邊界元法只是在解決線性問(wèn)題方面比較成熟，對(duì)于非線性問(wèn)題卻遠(yuǎn)非如此。如彈塑性力學(xué)問(wèn)題，雖然從四十多年前就有不少學(xué)者開(kāi)始對(duì)該課題進(jìn)行深入的研究［64，65］，但直到2000年，高效偉和Davies［66］才徹底解決了彈塑性邊界元法中強(qiáng)奇異域積分的計(jì)算問(wèn)題，并于2002年公開(kāi)發(fā)表了國(guó)際上第一個(gè)彈塑性邊界元程序［13］。我國(guó)的其他學(xué)者也對(duì)彈塑性邊界元法的發(fā)展做出了重要的貢獻(xiàn)［67��70］。
用邊界元法解決非線性和非均質(zhì)問(wèn)題時(shí)，由于很難求得控制方程的基本解，所以不得不用對(duì)應(yīng)于線性和均勻介質(zhì)問(wèn)題的基本解來(lái)建立積分方程，由此導(dǎo)致了域積分的出現(xiàn)。為了計(jì)算域積分，傳統(tǒng)的方法是將計(jì)算域離散成內(nèi)部網(wǎng)格，采用像有限元法中的方式計(jì)算域積分［64��70］，這樣就消除了邊界元法只需將邊界離散成單元的優(yōu)點(diǎn)。正是這種求解非線性和非均質(zhì)問(wèn)題中需要內(nèi)部網(wǎng)格的致命弱點(diǎn)，嚴(yán)重影響了邊界元法的發(fā)展。
為了避免使用內(nèi)部網(wǎng)格，國(guó)際上不少學(xué)者致力于無(wú)內(nèi)部網(wǎng)格邊界元法的研究，其中最常用的方法是將出現(xiàn)在積分方程中的域積分轉(zhuǎn)換成邊界積分。其中應(yīng)用最廣泛的方法是Breb