边界元法

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