《高等边界元法:理论与程序》共9章,第｡章为绪论,第2章介绍必要的数学知识,第3~6章介绍与位势问题相关的边界元法,第7~8章介绍线性和非线性力学问题的边界元法,第9章介绍求解多种介质问题的新方法｡《高等边界元法:理论与程序》展示了作者多年来的研究成果,如:将任意域积分转换成边界积分的径向积分法､求解大型非对称稀疏矩阵方程的同时消元回代法､计算弱奇异和近奇异积分的单元子分法､计算超奇异积分的幂级数展开法､建立一般非线性问题积分方程的源点隔离法以及用单一积分方程求解多种介质问题的界面积分方程法｡。

第1章绪论
1.1数值方法概述随着科学技术的不断发展，需要解决的工程问题也越来越复杂，对于大多数问题，由于求解问题几何形状的复杂性或计算介质的非线性，人们已经很难得到问题的解析答案。另外，随着计算机性能的日益提高，求解工程问题的数值计算方法也不断成熟，现在几乎所有的大型工程问题都要借助数值计算进行分析或评估，为科技人员的工程设计提供依据或参考。
目前，已经发展起来的用得较多的数值计算方法有五大类：有限差分法、有限体积法、有限单元法、无网格法和边界元法。
(1) 有限差分法(finite difference method)［1］。有限差分法是将所考虑的区域划分成网格，用差分近似微分，把微分方程变成差分方程。也就是通过数学上的近似，把求解微分方程的问题变换成求解关于节点未知量的代数方程问题。该方法简单、易懂，便于复杂微分方程的求解，因此在流体力学领域被广泛采用［2］。但当求解问题的几何形状复杂时，按空间坐标的离散变得困难，网格的正交性不易保留，导致计算精度降低。
(2) 有限体积法(finite volume method)［3］。有限体积法是基于物理问题控制方程的积分形式的数值方法。解题思路是：把计算域离散成有限个互不重叠的控制体单元（网格），通过将积分形式的控制方程作用于每一个体单元来建立离散的代数方程组。有限体积法既有有限差分法的特点，又有有限单元法的特点，方法简单、几何适应性好，是现代计算流体力学中占主导地位的数值方法［4］。其缺点是物理量的空间导数相关量（如通量）是由周围体单元中心处的值决定的，因此精度较低，特别是靠近边界的通量值更难计算准确。
(3) 有限单元法(finite element method)［5］。有限单元法是通过变分原理建立含权函数（形函数）的体积分形式方程的方法。解题思路是：首先将计算域离散成有限个互不重叠的有一定规则的节点组成的单元，然后对每个单元积分并通过组集形成总体代数方程组。有限单元法的单元可以按不同的连接方式进行组合，每个单元可以有不同的形状和材料性质，因此有限单元法具有几何适应性强、可灵活处理不同物性参数的优点，在各个领域都得到了广泛应用［6］。有限单元法的缺点是：物理量的空间导数是通过对形函数的求导得到的，因此精度要比物理量本身低一阶；在金属成形、优化计算、渗流问题自由面的确定等运动边界问题中，有限元网格可能产生畸变和重叠，以致计算精度下降或计算中止。
(4) 无网格法(meshless method)［7］。无网格法是基于构造点插值函数的数值方法，因此只需要在计算域内布置一系列的离散点即可，不需要网格单元，具有很强的解题灵活性［8，9］。但无网格法发展得还不够成熟，缺少坚实的理论基础和严格的数学证明，因此计算精度、守恒性等一直没有明确的答案。此外，无网格法是基于点的算法，因此布点数量和方案会影响计算时间和精度。另外，由于不使用单元，对于几何较复杂的运动边界问题，边界变化时要判断重新分布后的点是内部点、外部点还是边界点，这时会存在困难。
(5) 边界元法(boundary element method，BEM)［10］。边界元法是基于格林公式和问题的基本解将控制微分方程转化为边界积分方程的一种数值方法。其主要优点是：①只需要将边界离散成单元，因而准备数据简单、便于复杂几何问题的建模［11］；②能够自动满足无限远处的边界条件，因而适合于求解无限与半无限域问题［12］；③所求物理量的空间导数的计算公式可以解析地从基本边界积分方程中导出， 因此所求与导数相关的物理量（如通量、应力等）与物理量本身具有同样级别的精度［13］；④在求解运动边界问题时，移动边界节点的位移与原边界节点的坐标相加就自然形成了新的边界单元信息，不需要专门重构单元，也不会有网格重叠的问题［1416］。此外，由于在计算域的边界上有单元信息可用，所以通过单元积分很容易判断一点是内部、外部还是边界点。基于这些优点，边界元法在一些领域（如运动边界［1517］、裂纹［18，19］、接触［20，21］、辐射［22，23］、超薄结构［2428］、无限域和半无限域［14，29，30］、声学［31］等）的应用优于有限元法。
1.2边界元法的发展历史
边界元法已发展成为求解工程与科学问题的常用数值分析方法之一。它是一种在经典的积分方程基础上，吸收了有限元法的离散技术而发展起来的数值方法。边界元法通过采用一个满足无限或半无限域场方程的奇异函数——基本解作为权函数，将基于问题控制方程的域积分方程转化为边界积分方程，并将边界离散成边界单元来求解边界未知量的数值解。
边界元法的产生可追溯到19世纪，当时有人提出了积分等式和位势理论的数学问题，把线性偏微分方程的边值问题转化为边界积分方程求解。1905年，Fredholm将积分方程应用于求解弹性力学问题［32］。1953年，Muskhelishvili［33］和Kellogg［34］分别将积分方程法用于求解结构力学和位势问题。1965年，Mikhlin［35］解决了积分方程理论中的奇异性问题，为积分方程法在工程中的应用开辟了道路。
20世纪60年代，电子计算机的发展开创了数值求解积分方程的新时代，积分方程法作为数值计算方法开始应用于实际问题。1963年，Jaswon［36］采用间接边界积分方程法，成功地求解了位势问题和弹性力学问题。这种方法的主要思想是沿边界配置一种虚设的点源密度函数，先确定密度函数，再求边界上的未知物理量。其缺点是待求的点源分布函数是虚构的，不具有明确的物理意义，因此求解过程需要两步完成。但它具有场量方程和场量梯度方程相互独立的优点，因而易于组合求解各种复杂边界条件的边值问题。60年代，一些苏联学者对积分方程尤其是奇异积分方程的理论作了更为深入的研究，为进一步应用边界积分方程方法开辟了道路。与此同时，高速大型计算机的出现及其硬件的迅猛发展使离散求解积分方程成为可能。到了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］。书中用加权余量法推导出了边界积分方程，并指出加权余量法是最普遍的数值方法，如果以开尔文（Kelvin）解作为权函数，从加权余量法可导出弹性力学问题的边界积分方程，通过将边界离散成单元的方法可数值求解积分方程。至此，边界元法这一名称得到了国际公认。
自1978年第一届国际边界元法会议后，边界元法会议几乎每年在世界各地举办。世界各国已从基本理论与方法的研究向深广领域发展，大量论文和专著先后问世。在此时期，边界元法在数学分析理论和数值积分方法的研究、边界元法的完善和应用范围的拓宽以及边界元法应用软件等方面均得到飞速发展。边界元法的应用遍及固体力学、流体力学、波动学、传热学、电磁学等学科领域。
我国关于边界元法的研究大约开始于1978年，当时杜庆华在国内首先开创了工程中边界元法的研究，开始跟踪国际上这一领域的最新进展，其研究领域主要在固体力学方面［42］。与此同时，王泳嘉［43］、郑颖人等［44］开始了岩土工程中的边界元法研究，杨德全和赵忠生［45］在流体力学边界元法方面开了先河。值得一提的是，岑章志［46］对我国的弹塑性边界元法做了开拓性的研究工作，高效伟和Davies［13］发表了国际上第一个弹塑性边界元程序，结束了三十多年来在非线性力学边界元分析方面只有论文发表、没有程序公布的局面。
近年来，我国学者在快速多极边界元法研究方面取得了一系列的研究成果，代表性工作可见姚振汉与王海涛的著作［11］。此外，我国学者在近奇异积分计算［2428］、多种介质问题［4749］以及与CAD技术结合解决工程问题［50］等方面的研究工作也引人注目。边界元法的研究目前在我国正处于上升阶段［51］，近年来在一些回国学者的带领下越来越多的专家学者投身于边界元法的研究中。
1.3边界元法中的难点问题及其研究进展
边界元法在解决接触、断裂力学、运动边界、无限域与半无限域以及超薄结构等问题中具有独特的优势，被大量地用于科学与工程问题的计算分析，在许多应用领域，边界元法被公认为有限元法的一个重要补充。然而，传统的边界元法有下述几方面的弱点：
（1） 奇异性问题。边界元法中所用的基本解是奇异函数，数值计算时首先需要消除积分奇异性才能得到精确的计算结果。多年来已有大量的文献在计算效率与精度方面提出了不少计算奇异积分的方法，如线性单元的解析消除法［10，52］、弱奇异积分的单元子分法［13，53，54］、强奇异积分的间接计算法［13，55］、高阶奇异积分的直接计算法［5658］等。这些方法在计算弱奇异和强奇异积分时非常有效，但在处理高阶奇异积分时的稳定性还需要进一步考查。
（2） 满系数矩阵问题。边界元法在求解问题时所形成的方程组的系数矩阵是满阵，因而占有较大的计算机内存，难以解决大型工程问题。为了解决此问题，研究人员提出了两种有效的解决方法，一种是快速多极法，另一种是区域分解法。快速多极边界元法［11，59］，通过将奇异核函数进行级数展开的技术，降低矩阵向量相乘操作的计算量级和存储量级，达到节省内存和提高计算速度的目的。区域分解法［10，60］的基本思想是把总求解域划分成多个子域，对每个子域建立边界元矩阵方程，然后利用子域间公共节点上的面力平衡条件和位移相容性条件把子域方程组集成总体系统方程组。这样形成的系数矩阵是块状稀疏矩阵，可利用现有成熟的稀疏矩阵求解器（如LU分解法和GMRES迭代法（广义最小残量法））对系统方程组进行有效求解。区域分解法是边界元法中被广泛应用的方法，不仅能解决满系数矩阵问题，而且能通过按照材料性质划分子域的手段求解由不同材料组成的复合介质问题［61］，还能通过沿裂纹面划分子域的方法求解断裂力学问题［19］。最近，作者基于区域分解法，提出了求解非均匀介质问题的三步变量凝聚技术［49］，形成的系统方程组只有公共节点位移作为未知量，可以求解大型工程问题。虽然这些方程组的组建与求解技术能解决相当规模的工程问题，但对于超大型问题，系统方程组的求解仍然是一个极具挑战性的问题。最近，高效伟等提出了求解非对称稀疏方程组的同时消元回代法求解技术［62，63］，在计算效率和储存空间方面都有了显著的提高。
（3） 非线性和非均质问题中的域积分问题。传统的边界元法只是在解决线性问题方面比较成熟，对于非线性问题却远非如此。如弹塑性力学问题，虽然从四十多年前就有不少学者开始对该课题进行深入的研究［64，65］，但直到2000年，高效伟和Davies［66］才彻底解决了弹塑性边界元法中强奇异域积分的计算问题，并于2002年公开发表了国际上第一个弹塑性边界元程序［13］。我国的其他学者也对弹塑性边界元法的发展做出了重要的贡献［6770］。
用边界元法解决非线性和非均质问题时，由于很难求得控制方程的基本解，所以不得不用对应于线性和均匀介质问题的基本解来建立积分方程，由此导致了域积分的出现。为了计算域积分，传统的方法是将计算域离散成内部网格，采用像有限元法中的方式计算域积分［6470］，这样就消除了边界元法只需将边界离散成单元的优点。正是这种求解非线性和非均质问题中需要内部网格的致命弱点，严重影响了边界元法的发展。
为了避免使用内部网格，国际上不少学者致力于无内部网格边界元法的研究，其中最常用的方法是将出现在积分方程中的域积分转换成边界积分。其中应用最广泛的方法是Breb

前言
第1章绪论1
1.1数值方法概述1
1.2边界元法的发展历史2
1.3边界元法中的难点问题及其研究进展4
1.4本书重点内容概述7
参考文献8
第2章数学基础13
2.1符号标记13
2.1.1指标符号与求和约定13
2.1.2特殊的指标符号——δ14
2.1.3矢量简介15
2.2狄拉克δ函数15
2.2.1狄拉克δ函数的定义与性质16
2.2.2微分方程基本解17
2.3高斯公式与格林公式17
2.3.1高斯公式17
2.3.2格林公式18
2.4径向积分法18
2.5加权余量法22
2.5.1加权余量法一般公式22
2.5.2加权余量法算例演示24
2.6高斯数值积分公式27
2.6.1一维高斯积分28
2.6.2二维和三维高斯积分29
2.7大型线性方程组的求解方法——同时消元回代法29
2.7.1同时消元回代法求解线性方程组的基本公式30
2.7.2SEBSM求逆矩阵32
2.7.3基于SEBSM的线性方程组残差迭代求解法33
2.7.4同时消元回代法Fortran程序介绍35
2.7.5SEBSM应用示例39
2.7.6SEBSM与其他方法的比较42
附录2A选主元法求解稀疏矩阵方程的子程序SEBSM_SPARSE_PIVOT43
参考文献46
第3章位势问题47
3.1位势问题的积分方程47
3.2位势问题的基本解49
3.2.1二维位势问题基本解求解方法49
3.2.2三维位势问题基本解求解方法50
3.3位势问题的边界积分方程51
3.4位势梯度的积分方程54
3.5边界积分方程的离散54
3.5.1常单元55
3.5.2线性单元56
3.5.3高次单元61
3.6奇异积分和近奇异积分的计算65
3.6.1积分奇异性的定义66
3.6.2弱奇异积分的计算66
3.6.3强奇异积分的计算——常位势法70
3.6.4超奇异积分的计算——幂级数展开法71
3.6.5近奇异积分的计算——单元子分法79
7非齐次位势问题的积分方程法82
3.7.1泊松方程的边界-域积分方程82
3.7.2域积分到边界积分的转换方法83
3.8角点问题96
3.8.1混合单元法87
3.8.2多节点补充方程法87
3.8.3非连续单元法87
3.9程序介绍及算例90
3.9.1位势问题程序PTBEM介绍及算例90
3.9.2高阶奇异边界积分计算程序SIEPPEM介绍及算例96
附录3A二维常单元影响系数解析表达式102
附录3B面单元形函数106
附录3C三维面单元外法向向量n和雅可比行列式J(ξ，η)的计算式107
参考文献108
第4章热传导问题110
4.1瞬态常系数热传导问题111
4.1.1瞬态热传导问题的积分方程111
4.1.2域积分到边界积分的转换112
4.1.3积分方程的离散116
4.1.4用时间推进法求解瞬态热传导问题117
4.2稳态变系数热传导问题118
4.3瞬态变系数热传导问题120
4.3.1边界元基本方程120
4.3.2代数方程组的求解122
4.4稳态非线性热传导问题123
4.4.1非线性热传导问题的积分方程124
4.4.2域积分到边界积分转换的解析表达式124
4.4.3系统方程组的组集127
4.4.4系统方程组的迭代求解129
4.5瞬态非线性热传导问题129
4.5.1边界元基本方程129
4.5.2代数方程组的求解130
4.6内部热流积分方程131
4.6.1常系数问题内部热流积分方程131
4.6.2变系数问题内部热流积分方程132
4.6.3非线性问题内部热流计算公式136
4.7各向异性稳态热传导问题137
4.7.1各向异性热传导问题积分方程的一般关系式137
4.7.2基于格林函数基本解的各向异性积分方程关系式139
4.7.3各向异性材料温度梯度计算公式140
4.8程序介绍及算例143
4.8.1程序介绍143
4.8.2常系数瞬态问题数值算例149
4.8.3变系数问题数值算例152
4.8.4非线性问题数值算例158
参考文献165
第5章热辐射问题168
5.1热辐射简介170
5.1.1基本概念170
5.1.2黑体辐射173
5.1.3实际物体表面的辐射175
5.1.4介质内的辐射能量传递179
5.2表面辐射换热问题的边界元法181
5.2.1边界积分方程182
5.2.2边界积分方程的数值执行182
5.2.3奇异积分的处理183
5.3介质辐射换热问题的边界元法185
5.3.1边界-域积分方程185
5.3.2域积分到边界积分的转换187
5.3.3代数方程及求解188
5.3.4非灰体辐射换热问题189
5.4程序介绍及算例191
5.4.1程序介绍191
5.4.2算例分析194
参考文献202
第6章导热-辐射耦合换热问题204
6.1稳态线性导热-辐射耦合换热问题的边界元法205
6.1.1积分方程205
6.1.2积分方程的离散206
6.1.3导热-辐射耦合问题的系统方程组及其求解207
6.2瞬态非线性导热-辐射耦合换热问题的边界元法209
6.2.1边界积分方程209
6.2.2积分方程的离散210
6.2.3瞬态非线性方程组的求解212
6.3程序介绍及算例213
6.3.1程序介绍213
6.3.2算例分析215
参考文献225
第7章弹性力学问题227
7.1弹性力学基本方程227
7.1.1基本概念227
7.1.2基本方程229
7.1.3平面问题231
7.2弹性力学问题的基本解233
7.2.1位移基本解234
7.2.2应变基本解234
7.2.3应力基本解235
7.2.4面力基本解235
7.3弹性力学问题的边界积分方程236
7.3.1Betti互易功定理236
7.3.2Somigliana恒等式237
7.3.3边界积分方程238
7.4弹性力学问题边界积分方程的离散240
7.4.1边界离散及变量表示240
7.4.2代数方程组242
7.4.3边界积分的高斯公式计算243
7.5强奇异积分的处理244
7.5.1有限域问题244
7.5.2无限域问题245
7.5.3半无限域问题245
7.6边界应力计算247
7.6.1三维问题边界应力计算式247
7.6.2二维问题边界应力计算式249
7.7体积力的处理250
7.7.1重力251
7.7.2离心力252
7.8热弹性力学问题252
7.8.1基本方程252
7.8.2位移积分方程253
7.8.3内部应力积分方程253
7.8.4边界应力计算式256
7.8.5域积分到边界积分的转换257
7.9程序介绍及算例258
7.9.1程序介绍259
7.9.2算例分析262
附录7A弹性力学边界元程序BEMEL输入文件变量说明269
参考文献272
第8章非均质和非线性力学问题274
8.1一般固体力学问题的积分方程274
8.1.1一般固体力学问题积分方程的建立——源点隔离法274
8.1.2基于开尔文基本解的关系式276
8.1.3正则化的位移积分方程277
8.1.4应力积分方程278
8.2变系数各向同性弹性力学问题281
8.2.1位移积分方程282
8.2.2应力积分方程282
8.2.3域积分到边界积分的转换以及常用的径向积分解析计算式283
8.2.4变系数问题系统方程组290
8.3各向同性损伤弹性力学问题291
8.3.1损伤变量演化方程291
8.3.2系统方程组及牛顿-拉弗森迭代求解公式292
8.4各向异性弹性力学问题293
8.4.1各向异性弹性本构张量293
8.4.2弹性本构张量的坐标变换形式295
8.5弹塑性力学问题295
8.5.1弹塑性本构方程295
8.5.2新型弹塑性力学边界-域积分方程297
8.5.3弹塑性应力-应变本构张量空间导数的计算301
8.5.4弹塑性系统方程的迭代求解302
8.6程序介绍及算例305
8.6.1程序介绍305
8.6.2算例分析311
附录8A网格单元形函数319
8A.18节点三维网格单元形函数319
8A.220节点三维网格单元形函数320
参考文献321
第9章多种介质问题322
9.1热传导问题多种介质界面积分方程322
9.1.1基于“域积分界面退化”技术的界面积分方程323
9.1.2分片各向同性介质界面积分方程327
9.2固体力学问题多种介质界面积分方程328
9.2.1基于“域积分界面退化”技术的界面积分方程328
9.2.2分片各向同性介质弹性力学界面积分方程332
9.3基于“三步变量凝聚法”的多域边界元法334
9.3.1子域划分与节点分类335
9.3.2固体力学问题三步变量凝聚法335
9.3.3瞬态热传导问题三步变量凝聚法338
9.4程序介绍及算例342
9.4.1程序介绍343
9.4.2多种介质热传导问题算例346
9.4.3多种介质弹性力学问题算例350
附录9A自由项系数fijs的正则化计算公式356
附录9B自由项系数fijmns的正则化计算公式359
参考文献360