第8章 边界元法

件，但不能精确地满足微分方程。
现在从一般性的加权余量法展开讨论，假设定解问题为式（8-9a）、式
（8-9b）和式（8-9c）所描述的三维线性泊松场。设其近似解 是某一线性
无关的完备函数集合
在一般情况下，把近似解 代入该定解问题，微分方程（8-9a）和边界条件 （8-9b）、（8-9c）都将不能精确满足，由此产生相应的误差，其余量可分 别表示为
上式称为格林第一公式。若将ψ和φ交换位置，即对向量φ∇ψ进行同样的处 理，便得
第8章 边界元法
以式（8-3）减去式（8-4），则有
上式称为格林第二公式，亦称为格林定理。
8.2.2 基本解
若考虑一线性微分方程
式中，L是线性微分算子，f是给定的激励源。则满足方程
的解u（r，r′）称为对应于方程（8-6）的基本解。式（8-7）中的激励源项为 狄拉克δ函数，由定义式（7-20）可见其具有点源性质。u（r，r′）亦可称为 下列方程的基本解，即
已如前述，从数学意义上分析，加权余量法是其他多种数值计算方法的 基础，取决于不同的权函数W的选择，可派生出不同类型的相应计算方法。 就边界元法而言，即可直接由加权余量法出发，导得构造边界元法的数学基 础———边界积分方程，并选取相应的权函数为基本解展开阐述。
第8章 边界元法
8.3 边界积分方程
8.3.1 边界积分方程
基于式（8-17），还可导出应用于边界元法的间接边界积分方程（间接 法公式）［8］。但就边界元法而言，直接法比间接法的计算步骤少，计算 精度高，故本书以直接边界积分方程为基础，展开叙述边界元法。
第8章 边界元法
8.4 边界元方程及方法实施
在给定边界条件和场域几何形态的情况下，采用解析的方法求解边界积 分方程是十分困难的，因此，作为一种有效的数值计算方法———边界元法， 借助于有限元技术，通常可由以下步骤组成：
第8章 边界元法
如前所述，为了使这些在场域内和边界S1、S2上的余量为最小，可引入一个 权函数W，使之在平均意义上令余量的加权积分为零。根据误差分布原理 ［5］ ，不难导得
上式表明，所选择的近似解既不满足基本方程，也不满足相应的边界条件。 因此，式（8-14）可以看作是前述加权余量式（7-4）的推广，并由此可以求 出近似解 。
2）方程组阶数降低，输入数据量减少。如前所述，待求量将仅限于边界节 点，这不仅简化了问题的前处理过程，而且大幅度降低了待求离散方程组 的阶数。
3）计算精度高。本方法直接求解的是边界广义场源的分布。根据不同的问 题，广义场源可以是位势、场源或等效场源。场域中任一点的场量将通过 线性叠加各离散的广义场源的作用而求得，毋需再经微分运算。此外，由 于只对边界离散，离散化误差仅仅来源于边界。所以边界元法较之有限元 法，可望有较高的计算精度。
显然，如静电场中泊松方程的基本解［式（8-11a）］，即表示在无界空 间位矢为r′的点上放置一电量为ε0的正电荷，它在与其相距r处所产生的电位 值φ=1/(4πr)。由此可知，呈体电荷密度ρ分布的场源在该场点产生的电位就 等于此基本解乘以ρdV′/ε0，然后对应于源区的体积分，即
第8章 边界元法 8.2.3 加权余量法的推广
式（8-34）又可改写成
第8章 边界元法
因为在边界L1上有N1个单元属于第一类边界条件，即其N1个单元上的u 值是已知的，但其q值未知；而边界L2上对应的N2（=N-N1）个单元属于第二 类边界条件，即其N2个单元上的q值已知，但u值未知。因此，离散的边界积 分方程的未知量应由N1个q值和N2个u值所组成。式（8-36）是对应于第i个节 点所列出的离散边界积分方程，就整体N个边界节点的集合而言，即构成N 阶方程，可写成如下矩阵形式：
为了书写简便起见，往后将近似解 改记为u。从而，式（8-14）可以重写为
根据格林第二公式［式（8-5）］，Βιβλιοθήκη Baidu式左边可表示为
式中，在边界S1、S2上，记 15），经整理得
。将式（8-16）代入式（8-
上式是电磁场边界积分方程的原始公式，由此可推导出直接边界积分方程。
第8章 边界元法 8.3.2 直接边界积分方程（直接法公式）
本章最后给出二维边界元法典型应用的示例，讨论了方法实施的全过 程，以及可供参考使用的计算程序，并通过本方法与其他方法计算结果的 对比，进一步展示了本方法的特点。
8.1 概述
边界元法（Boundary Element Method，简称BEM）是近20余年来发展 形成的一种数值计算方法。该方法的工程应用起始于弹性力学，现进而应用 于流体力学、热力学、电磁工程、土木工程等诸多领域，并已从线性、静态 问题延拓到非线性、时变问题的研究范畴。
第8章 边界元法 8.4.1 常数单元
常数单元是指每个边界单元上的u和q值都设定为相应的常数，且等于该 单元中点上的值。各单元中心即其两端点连线的中心点，亦称节点，如图83所示。图中L1、L2分别标记给定的第一类和第二类边界条件所对应的边界。
设场域D内位函数u满足拉普拉斯方程，则直接边界积分方程（8-28）可以写为
将式（8-18）代入式（8-17），该式左边
式中，ui是V域内节点i处的u值。因此，式（8-17）可以写成
由上式可见，一旦求出边界上的物理量u 点的物理量值ui。
，便可解得V域内任一
第8章 边界元法
可以看出，当求解边界上的物理量时，在场点 与源点重合（即r=0）处，式（8-19）中的面积分 项会出现奇异积分。此时处理方法如下：
5）同样基于边界积分方程，在上述边界元法所得离散解的基础上，可得场 域内任一点的位函数与场量解。
本节讨论应用于二维问题的边界元法。关于三维问题的边界元法，其 基本思想类同，但由于离散的边界单元将是平面或曲面形单元，处理过程 较为繁复，限于篇幅，不再展开阐述和讨论。
第8章 边界元法
二维场的边界积分方程已由式（8-28）给出。该二维场域D的边界L是 一维曲线，现按有限元离散方法，将边界离散成N个边界单元（L1，L2，…， LN），并规定单元序号（或节点序号 ）与边界定向线段L的走向一致，即所 论场域D始终位于L的左侧。如图8-2所示。插值函数有各种类型，基本上可 分为常数型、线性型和高次插值。下面从最简单的常数单元入手，推导边界 元方程。
重新排列上式，将所有包含有未知量的项移置方程的左端，而将已知项置于方 程的右端，可得重排后的N阶线性方程组，即边界元方程为
式中，X表示由未知量u和q所组成的列向量；F是N维列向量，表示给定的边 界条件； A为N×N阶系数矩阵，表征了节点i与各单元j之间的关联。一旦方程 （8-38）解出，即可求得边界上所有未知的u和q值，而按式（8-29）场域内任 一点的位函数u的计算公式为
设边界面S1光滑，在该边界面上，以场点i为 球心，半径r0=ε作半球面 ，如图8-1所示。然 后令ε→0，以求得相应面积分在点i上的极限值。 这一分析将包含以下三种情况： （1）场点位于S1面外，且场点不在V域内 此时，基本解满足拉普拉斯方程∇2W=0，式（8-17）可简化成
（2）场点在S1面上
第8章 边界元法
限于篇幅，这里不讨论基本解的求解方法，而是直接给出电磁场工程问题， 如静态电磁场问题常用的基本解。
第8章 边界元法
静态场问题可由泊松方程或拉普拉斯方程的定解问题一般地描述为
其二维问题的基本解为 三维问题的基本解为
第8章 边界元法
式中，r是源点到场点间的距离；u则代表位势或场量的某一分量。
从以上基本解的定义可以看出，基本解的实质是集中量（点源） Cδ（r-r′）在空间产生的效应。就线性微分方程而言，如果激励场源是一连续 分布量，那么它所产生的效应可以根据线性叠加原理，表示成无数个集中量 所产生的效应的叠加。也就是说，连续分布量所产生的效应可以用基本解乘 以连续分布量的密度函数的积分来表示。
第8章 边界元法
8.2 基础知识
8.2.1 格林公式
设V为空间中某一闭域，其表面为S。若有两个标量函数φ和ψ，它们在V
域内及S面上分别存在连续的一阶和二阶偏导数，则所构成的向量ψ∇φ 满
足如下的高斯散度定理：
式中，en为S面的外法线方向的单位向量； 为法向导数。根据向量恒等式
将式（8-2）代入式（8-1）可得
在第7章中已经讨论了可以构成矩量法、伽辽金有限元法等的共同数学
基础———加权余量法。该方法表明，给定微分方程的近似解 在场域内不
能精确地满足微分方程，因而存在余量
，于是通过令该余量在平均意
义上，其加权积分为零，即得加权余量式（7-4）。应该指出，该式对应的
是加权余量法的最简情况，即所选择的近似函数 可以精确地满足边界条
直接边界积分方程中的未知量是边界上客观的物理量。如位函数φ、A， 磁场强度H及电场强度E等。一旦这些未知量被确定，场域内任一点上的物理 量值便即可求得。这就是应用于边界元法的直接法。 在电磁场问题中，现取权函数W为基本解。仍以三维泊松场为例，由式 （8-11a）可知，基本解W=1/(4πr)满足以下方程：
当ε→0时，上式右边第一项中
，故有
且以基本解代入式（8-21）右边第二项，并注意到当ε→0时，
式中，Ω为场点 i 对于 面所张的立体角 将式（8-22）、式（8-23）代入式（8-21），可得
，得
第8章 边界元法
（3）场点在S2面上 类同于场点在S1面上时的分析，可以导得
显然，以上三种情况可以统一表示为 式中
第8章 边界元法
式（8-26）亦可表示为
若为拉普拉斯方程定解问题，则f=0， 对于二维场问题，通过类似的推导，最终统一表达式应为
式中
且其基本解
。
第8章 边界元法
式（8-27）和式（8-28）即为直接边界积分方程。当边界面（线）光滑， 且场点i位于边界上时，对应于三维和二维问题的ci=1/2（Ω=2π；θ=π）。由 此，可求出边界上的未知量。然后，再令直接边界积分方程中的ci=1，即可 解出场域内任一点处的场量。
1）边界S被离散成一系列边界单元，在每个单元上，假定位势及其导数是 按节点值的内插函数形式变化。
2）基于边界积分方程，按边界单元上节点的配置，在相应节点上建立离散 方程。
3）采用数值积分法，计算每个单元上的相应积分项。
4）按给定的边界条件，确立一组线性代数方程组，即边界元方程。然后， 采用适当的代数解法，解出边界上待求的位势或其导数的离散解。
4）易于处理开域问题。本方法只对有限场域或无限场域的有限边界进行离 散化处理并求解，因此特别适用于开域问题。
第8章 边界元法
然而，边界元法与有限元法相比较，其明显的不足之处是： 1）系数矩阵为非对称性的满阵。显然，这就引发了应用计算机求解大型离散 方程组的困难，从而约束了边界元方程组的阶数。 2）系数矩阵元素值需经数值积分处理，故系数矩阵的建立需要较多的计算机 时。 3）不易处理多种媒质共存的问题。
第8章 边界元法
边界元法是把边值问题等价地转化为边界积分方程问题，然后利用有限 元离散技术所构造的一种方法，其主要特点是：
1）降低问题求解的空间维数。本方法将给定场域的边值问题通过包围该场 域边界面上的边界积分方程来表示，从而降低了问题求解的空间维数。也 就是说，三维问题可利用边界表面积分降维为二维问题；而二维问题则利 用边界的线积分降维为一维问题。因此，有限元离散仅对应于二维曲面单 元或一维曲线单元，使方法的构造大为简化。
第8章 边界元法
第8章 边界元法
本章基于加权余量法，阐述了构造边界元法的数学基础———边界积 分方程。由此，引入微分方程基本解和格林公式，进一步导出了对应于边 界上未知量为场量φ、A、E、和H的直接边界积分方程。
以数值求解边界积分方程为目的，本章介绍了两种最基本的边界元法的 构造模式：常数单元和线性单元的计算模式。并借助于高斯求积公式给出各 系数矩阵的元素值。
第8章 边界元法
当边界离散后，按边界单元上节点的配置，上式可改写为
式中，i为节点序号；j为单元序号。由于在各个边界单元Lj（j=1，2，…，N） 上u与q均分别设定为相应的常数，故可将其提出积分号，得 各单元Lj上的积分仅与节点i和单元j相关。令
和
第8章 边界元法
和Gij一般可由数值积分算出，对于边界几何形状非常简单的情况，当 然也可以有解析解。这样，式（8-30）即为 前已指出，惟有当场点与源点重合时，即i=j时，ci=1/2（边界光滑时），其 余均为零。故若再令