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