第8章边界元法ppt课件.ppt

《第8章边界元法ppt课件.ppt》由会员分享，可在线阅读，更多相关《第8章边界元法ppt课件.ppt（48页珍藏版）》请在三一办公上搜索。

1、第 8 章 边 界 元 法,本章基于加权余量法，阐述了构造边界元法的数学基础边界积分方程。由此，引入微分方程基本解和格林公式，进一步导出了对应于边界上未知量为场量、A、E、和H的直接边界积分方程。,以数值求解边界积分方程为目的，本章介绍了两种最基本的边界元法的构造模式：常数单元和线性单元的计算模式。并借助于高斯求积公式给出各系数矩阵的元素值。,本章最后给出二维边界元法典型应用的示例，讨论了方法实施的全过程，以及可供参考使用的计算程序，并通过本方法与其他方法计算结果的对比，进一步展示了本方法的特点。,8.1 概述,边界元法（Boundary Element Method，简称BEM）是近20余年

2、来发展形成的一种数值计算方法。该方法的工程应用起始于弹性力学，现进而应用于流体力学、热力学、电磁工程、土木工程等诸多领域，并已从线性、静态问题延拓到非线性、时变问题的研究范畴。,边界元法是把边值问题等价地转化为边界积分方程问题，然后利用有限元离散技术所构造的一种方法，其主要特点是：,1）降低问题求解的空间维数。本方法将给定场域的边值问题通过包围该场域边界面上的边界积分方程来表示，从而降低了问题求解的空间维数。也就是说，三维问题可利用边界表面积分降维为二维问题；而二维问题则利用边界的线积分降维为一维问题。因此，有限元离散仅对应于二维曲面单元或一维曲线单元，使方法的构造大为简化。2）方程组阶数降低

3、，输入数据量减少。如前所述，待求量将仅限于边界节点，这不仅简化了问题的前处理过程，而且大幅度降低了待求离散方程组的阶数。3）计算精度高。本方法直接求解的是边界广义场源的分布。根据不同的问题，广义场源可以是位势、场源或等效场源。场域中任一点的场量将通过线性叠加各离散的广义场源的作用而求得，毋需再经微分运算。此外，由于只对边界离散，离散化误差仅仅来源于边界。所以边界元法较之有限元法，可望有较高的计算精度。4）易于处理开域问题。本方法只对有限场域或无限场域的有限边界进行离散化处理并求解，因此特别适用于开域问题。,然而，边界元法与有限元法相比较，其明显的不足之处是：,1）系数矩阵为非对称性的满阵。显然

4、，这就引发了应用计算机求解大型离散方程组的困难，从而约束了边界元方程组的阶数。2）系数矩阵元素值需经数值积分处理，故系数矩阵的建立需要较多的计算机时。3）不易处理多种媒质共存的问题。,8.2 基础知识,8.2.1 格林公式,设V为空间中某一闭域，其表面为S。若有两个标量函数和，它们在V域内及S面上分别存在连续的一阶和二阶偏导数，则所构成的向量满足如下的高斯散度定理：,式中，en为S面的外法线方向的单位向量；为法向导数。根据向量恒等式,将式（8-2）代入式（8-1）可得,上式称为格林第一公式。若将和交换位置，即对向量进行同样的处理，便得,以式（8-3）减去式（8-4），则有,上式称为格林第二公式

5、，亦称为格林定理。,8.2.2 基本解,若考虑一线性微分方程,式中，L是线性微分算子，f是给定的激励源。则满足方程,的解u（r，r）称为对应于方程（8-6）的基本解。式（8-7）中的激励源项为狄拉克函数，由定义式（7-20）可见其具有点源性质。u（r，r）亦可称为下列方程的基本解，即,限于篇幅，这里不讨论基本解的求解方法，而是直接给出电磁场工程问题，如静态电磁场问题常用的基本解。,静态场问题可由泊松方程或拉普拉斯方程的定解问题一般地描述为,其二维问题的基本解为,三维问题的基本解为,式中，r是源点到场点间的距离；u则代表位势或场量的某一分量。,从以上基本解的定义可以看出，基本解的实质是集中量（点

6、源）C（r-r）在空间产生的效应。就线性微分方程而言，如果激励场源是一连续分布量，那么它所产生的效应可以根据线性叠加原理，表示成无数个集中量所产生的效应的叠加。也就是说，连续分布量所产生的效应可以用基本解乘以连续分布量的密度函数的积分来表示。,显然，如静电场中泊松方程的基本解式（8-11a），即表示在无界空间位矢为r的点上放置一电量为0的正电荷，它在与其相距r处所产生的电位值=1/(4r)。由此可知，呈体电荷密度分布的场源在该场点产生的电位就等于此基本解乘以dV/0，然后对应于源区的体积分，即,8.2.3 加权余量法的推广,在第7章中已经讨论了可以构成矩量法、伽辽金有限元法等的共同数学基础加权

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.3 边界积分方程,8.3.1 边界积分方程,为了书写简便起

9、见，往后将近似解 改记为u。从而，式（8-14）可以重写为,根据格林第二公式式（8-5），上式左边可表示为,式中，在边界S1、S2上，记。将式（8-16）代入式（8-15），经整理得,上式是电磁场边界积分方程的原始公式，由此可推导出直接边界积分方程。,8.3.2 直接边界积分方程（直接法公式）,直接边界积分方程中的未知量是边界上客观的物理量。如位函数、A，磁场强度H及电场强度E等。一旦这些未知量被确定，场域内任一点上的物理量值便即可求得。这就是应用于边界元法的直接法。,在电磁场问题中，现取权函数W为基本解。仍以三维泊松场为例，由式（8-11a）可知，基本解W=1/(4r)满足以下方程：,将式（

10、8-18）代入式（8-17），该式左边,式中，ui是V域内节点i处的u值。因此，式（8-17）可以写成,由上式可见，一旦求出边界上的物理量u和，便可解得V域内任一点的物理量值ui。,可以看出，当求解边界上的物理量时，在场点与源点重合（即r=0）处，式（8-19）中的面积分项会出现奇异积分。此时处理方法如下：,设边界面S1光滑，在该边界面上，以场点i为球心，半径r0=作半球面，如图8-1所示。然后令0，以求得相应面积分在点i上的极限值。这一分析将包含以下三种情况：,（1）场点位于S1面外，且场点不在V域内,此时，基本解满足拉普拉斯方程2W=0，式（8-17）可简化成,（2）场点在S1面上,当0时

11、，上式右边第一项中，故有,且以基本解代入式（8-21）右边第二项，并注意到当0时，得,式中，为场点 i 对于 面所张的立体角,将式（8-22）、式（8-23）代入式（8-21），可得,（3）场点在S2面上,类同于场点在S1面上时的分析，可以导得,显然，以上三种情况可以统一表示为,式中,式（8-26）亦可表示为,若为拉普拉斯方程定解问题，则f=0，,对于二维场问题，通过类似的推导，最终统一表达式应为,式中,且其基本解。,式（8-27）和式（8-28）即为直接边界积分方程。当边界面（线）光滑，且场点i位于边界上时，对应于三维和二维问题的ci=1/2（=2；=）。由此，可求出边界上的未知量。然后，再

12、令直接边界积分方程中的ci=1，即可解出场域内任一点处的场量。,基于式（8-17），还可导出应用于边界元法的间接边界积分方程（间接法公式）8。但就边界元法而言，直接法比间接法的计算步骤少，计算精度高，故本书以直接边界积分方程为基础，展开叙述边界元法。,8.4 边界元方程及方法实施,在给定边界条件和场域几何形态的情况下，采用解析的方法求解边界积分方程是十分困难的，因此，作为一种有效的数值计算方法边界元法，借助于有限元技术，通常可由以下步骤组成：,1）边界S被离散成一系列边界单元，在每个单元上，假定位势及其导数是按节点值的内插函数形式变化。2）基于边界积分方程，按边界单元上节点的配置，在相应节点上

13、建立离散方程。3）采用数值积分法，计算每个单元上的相应积分项。4）按给定的边界条件，确立一组线性代数方程组，即边界元方程。然后，采用适当的代数解法，解出边界上待求的位势或其导数的离散解。5）同样基于边界积分方程，在上述边界元法所得离散解的基础上，可得场域内任一点的位函数与场量解。,本节讨论应用于二维问题的边界元法。关于三维问题的边界元法，其基本思想类同，但由于离散的边界单元将是平面或曲面形单元，处理过程较为繁复，限于篇幅，不再展开阐述和讨论。,二维场的边界积分方程已由式（8-28）给出。该二维场域D的边界L是一维曲线，现按有限元离散方法，将边界离散成N个边界单元（L1，L2，LN），并规定单元

14、序号（或节点序号）与边界定向线段L的走向一致，即所论场域D始终位于L的左侧。如图8-2所示。插值函数有各种类型，基本上可分为常数型、线性型和高次插值。下面从最简单的常数单元入手，推导边界元方程。,8.4.1 常数单元,常数单元是指每个边界单元上的u和q值都设定为相应的常数，且等于该单元中点上的值。各单元中心即其两端点连线的中心点，亦称节点，如图8-3所示。图中L1、L2分别标记给定的第一类和第二类边界条件所对应的边界。,设场域D内位函数u满足拉普拉斯方程，则直接边界积分方程（8-28）可以写为,当边界离散后，按边界单元上节点的配置，上式可改写为,式中，i为节点序号；j为单元序号。由于在各个边界

15、单元Lj（j=1，2，N）上u与q均分别设定为相应的常数，故可将其提出积分号，得,各单元Lj上的积分仅与节点i和单元j相关。令,和,和Gij一般可由数值积分算出，对于边界几何形状非常简单的情况，当然也可以有解析解。这样，式（8-30）即为,前已指出，惟有当场点与源点重合时，即i=j时，ci=1/2（边界光滑时），其余均为零。故若再令,式（8-34）又可改写成,因为在边界L1上有N1个单元属于第一类边界条件，即其N1个单元上的u值是已知的，但其q值未知；而边界L2上对应的N2（=N-N1）个单元属于第二类边界条件，即其N2个单元上的q值已知，但u值未知。因此，离散的边界积分方程的未知量应由N1个

16、q值和N2个u值所组成。式（8-36）是对应于第i个节点所列出的离散边界积分方程，就整体N个边界节点的集合而言，即构成N阶方程，可写成如下矩阵形式：,重新排列上式，将所有包含有未知量的项移置方程的左端，而将已知项置于方程的右端，可得重排后的N阶线性方程组，即边界元方程为,式中，X表示由未知量u和q所组成的列向量；F是N维列向量，表示给定的边界条件；A为NN阶系数矩阵，表征了节点i与各单元j之间的关联。一旦方程（8-38）解出，即可求得边界上所有未知的u和q值，而按式（8-29）场域内任一点的位函数u的计算公式为,基于同样的离散化过程，其离散形式是,值得注意的是，与式（8-35）不同，现节点i位

17、于场域内部，不会出现i=j的情况，故ci=1。,若继续求解场域内点i处的场强，即u的导数时，由于式（8-39）中的被积函数只有基本解W与点i相关，即只有W是r的函数，所以,式中，=x，y。可以看出，与u有相同的精度，即场强与位势有同阶的计算精度，这是边界元法的固有特点。,若问题满足的是泊松方程，则有,将场域D离散化为M个面单元Dk（k=1，2，M），令,通常f为已知的场源分布项，因此Bi项的引入，仅使式（8-38）中右端项F有所变化，但并不增加未知量。式（8-42）的离散形式为,8.4.2 线性单元,线性单元是将每个边界单元的端点取为节点，如图8-4所示，边界L离散为M个线性单元（设定每一个单

18、元内u和q的数值呈线性变化）。,由于单元中u和q值不是常数，所以式（8-30）中的u和 不可能提到积分号外，此时系数矩阵的建立较常数单元费时，其方法如下：,取任意j号单元（Lj为单元长度），建立如图8-5所示的局部坐标系。此时单元上任一点的u和q值可用相关节点值uj、uj+1和qj、qj+1以及N1、N2两个线性函数来表达，即,式中，N1、N2为插值基函数,亦称为形状函数，其特点是：在该单元节点上相应取值为1；在其余单元的节点上取值为零。,将式（8-44）代入式（8-30），可得,令,则有,上式亦可表示为,式中,因为边界L是一闭合曲线，所以当j=1时，（j-1）即为M。式（8-50）还可简化为

19、,式中,对应于i=1，2，M所有节点列出的离散方程（8-52），可写成矩阵形式为,显然，上式与式（8-37）的形式完全相同，只是系数矩阵元素计算关系式不同。经整理后，同样最终可得如式（8-38）所示的边界元方程。,对于二次或更高次的边界单元来说，差异仅在于插值基函数，即形状函数的构造将更为复杂，这时对应的是曲线形元素（常数及线性单元属于直线单元）。原则上，它们对曲线边界的拟合将更好，计算精度高，但由于形状函数的复杂性所带来的数值积分误差较大，其结果往往得不偿失。因此，当边界几何形状较简单或剖分足够精细时，通常采用线性单元即已可满足分析需要。,8.4.3 系数矩阵元素的确定,在建立系数矩阵时，会

20、遇到诸如式（8-32）、式（8-33）和式（8-47）、式（8-48）等所示的积分计算。这里分别就常数单元与线性单元讨论如下：,（1）常数单元,1）主对角元素,根据式（8-32）,式中，被积函数是基本解W的梯度在单元i的法线方向en上的投影。可以看出，由于en与单位向量（节点i到单元i上任一点的单位向量，即W的梯度方向）相互垂直，故积分为零，即,将基本解式（8-10a）代入式（8-33），有,由于常数单元的节点，即等效源的源点，位于单元中心，因此，在0lLi/2区域内，l0与r0的方向相反，如图8-6所示。所以,2）非对角线元素,前面已经指出，一般情况下，式（8-32）和式（8-33）不可能应

21、用解析方法求积。因此借助于3.3节中一维高斯求积公式（3-11），可以方便地求得Hij和Gij如下：,和,式中，en是单元j的法线方向；角是en与x轴的夹角，如图8-7所示。可以表示为，Ak是权系数，N为高斯积分点数，k是第k个高斯积分点的坐标。,（2）线性单元,1）主对角线元素,采用线性单元离散曲线边界，相当于用N条直线段去逼近曲线。设定节点为线段的端点，显然节点处曲线不再光滑，如图8-8所示。根据式（8-51）和式（8-52），对角线元素为,和,式中，ci=1-/2，是节点i所张的平面角。由于单元中r0与en相互垂直，经过与式（8-54）和式（8-56）相类似的推导，可以得出,和,2）非对

22、角线元素,此时,根据式（8-47）和式（8-48），并注意到l=Lj（1+）/2，所以有,同理可得,和,式中，r是的函数，,如上所述，应用高斯求积公式，即可算得以上各系数矩阵的元素值。,8.5 典型算例,仍以长直接地金属槽中电位和电场强度分布的计算为例，具体阐述二维边界元法的实施及其计算精度的评价。该金属槽如图8-9所示，侧壁和底面电位为零，顶盖电位的相对值为100。基于场分布的对称性，此二维电场的求解区域可进而压缩而归结为D/2。这样，对应于本问题的数学模型为,采用8.4节中线性单元的离散化模式，沿计算场域边界按逆时针方向编写单元的节点号，如图8-10所示。然后，基于线性插值函数，建立系数矩

23、阵H和G，由所得边界元方程（8-38）解出边界上的未知量后，即可计算场域内任一点的电位与场强。以上采用线性单元边界元法的数值求解过程，可应用列于附录8.1的线性单元边界元法的通用计算程序（BEM2D）来实现。为此，按该程序输入原始给定数据的子程序（INPUT）的规定顺序与格式，依次输入以下数据：,在场域边界剖分完全相同的前提下，表8-1和表8-2给出当a=20，b=17时，分别应用有限元法、边界元法和解析法计算典型场点处电位及电场强度的分布，并列出所得数值解的相对误差。由表8-1、表8-2计算结果可以看出，边界元法和有限元法对于位函数的计算，都有令人满意的计算程度。但是，在计算场强时，很明显，边界元法较之于有限元法有更高的计算精度，这是因为采用边界元法计算时，如8.4.1节所述，场强与位函数有同阶的计算精度，而用有限元法求场强时，则需再进行一次微分运算，从而导致计算精度的降低。,