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