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