边界单元法基础(直接法)

1、边界单元法基础（直接法）一、概述近年来在边界法方面人们发表了大量的文章和著作。这些方法是以不同的名称而提出来的，如“边界积分方程方法”“边界积分解”，等等。这种方法的数值解形式是把所考虑的域的边界划分为一系列的单元。边界单元法简称BEM是七十年代兴起的一种新的计算方法。它将边界上的广义位移和广义力作为独立变量且同时用满足场方程的奇异函数（源函数）作为加权函数。所以，它是一种特殊格式的加权余量法。边界元法只需将求解域的边界划分成单元，故使求解问题的维数降低，如三维问题可转变成二维问题求解。二维问题可化为一维问题。因而，输入数据大为减少，计算时间缩短。由于它只对边界离散，故离散误差仅为来源于边界，

2、而域内变量可由解析式的离散形式直接求得。因此，提高了计算精度。求域内变量时，只须改变其数量和坐标位置即可。和有限元法一样，边界元法可广泛地用来解决各种工程问题，如弹性力学、断裂力学、塑性力学、流体力学、温度场和电磁场等。边界元法分为直接法和间接法。直接法是用物理意义明确的变量来建立积分方程，其中未知函数就是所求的物理量在边界上的值；间接法是用物理意义不一定很明确的变量来建立积分方程，如位势问题中用单层位势和双层位势表示物理量。本部分着重叙述直接法。在用加权余量法建立积分方程时，所使用的权函数是数学上的“基本解”。基本解在数学上是作为微分方程的特殊的非齐次解定义的，它在每个问题上分别具有不同的物

3、理含义。求这个解，特别是便于解析的形式，一般是不容易的，这是数学上的难点。然而，除了特殊问题以外，主要微分方程的基本解，数学教科书中有所推导，工程技术人员可直接引用。边界元法另一个问题是，代数方程组的系数矩阵一般是非对称的，且非零系数矩阵为满秩矩阵，这是由于边界点与全部边界单元有关得出的，编程序时需要注意这一点。我们先介绍位势问题的边界单元法公式。这些基本概念对任何其它工程问题是类似的。然后再介绍利用边界元法求解弹性体受力分析问题。这里强调的是方法在工程中的应用。作为有限元法数值法的一个补充。二、泊松方程的边界单元法1积分方程的建立和基本解为了说明边界单元法的积分方程是如何由加权余量法推导得来

4、的，我们以泊松方程为例来阐述其全部求解过程，这对了解其它问题的边界元法求解是有益的。考虑势函数。它在域内满足微分方程，即 （在上） （9-1）在边界上满足边界条件，即 （在上） （在上） （9-2）如图9-1所示，求解域的总边界。图9-1 位势问题的域和边界可以证明，对于泊松方程或拉普拉斯方程（），一般加权余量表示式为 （9-3）式中，W是权函数，在边界法中可令W，称为相应于方程（9-1）的基本解，后面描述它。对式（9-3）左边第二项进行分部积分，即 （9-4）式中右边第一项的被积函数形式称爱因斯坦求和约定，即把式（9-4）右端第一项再次分部积分，得所以把上式代入式（9-3）或去（），考虑到，

5、把左右两边界积分合并时，得 （9-5）这就是求解泊松方程的积分方程。如果f = 0，即为求解拉普拉斯方程的积分方程。下面来导出相应该方程的基本解。假如有一单位势作用在物体上的P点，那么基本解就是满足下列方程2*（p,q）+（p-q）=0 （9-6）的二点函数，其中p，q为无限域中任意两点，如图9-2所示。图9-2 无限域中的p和q点取p点为坐标原点，并把式（9-6）写成极坐标表达式，即（9-7）式中，是狄拉克函数，它具有如下性质： （9-8）r为p点到q点之间的距离。方程（9-6）或（9-7）的解称为拉普拉斯方程的基本解。在式（9-7）中，当时，则为 （9-9）它的解为 （9-10）式中，C和

6、2为积分常数，可由r = 0时的条件确定。由于 （9-11）应用格林公式 （9-12）由图9-2（b）所示，以p为圆心，以为半径作一小圆，将式（9-10）代入式（9-12）并使时取极限，于是得由式（9-11）得 （9-13）由于是势函数，C1可取作零，所以 （9-14）式（9-14）即为二维拉普拉斯问题的基本解。同理可求得三维拉普拉斯问题的基本解为 （9-15）现在根据式（9-5）来建立积分方程。由于于是得 （9-16）式中；考虑时，可写成式（9-16）就是域内任一点p的函数值与边界的积分关系，除了微分方程非齐次项f的积分项外，其余都是边界积分表达式，故称积分方程。图9-3 用半球包围的边界点

7、下面我们来建立边界解的积分表达式。如果在式（9-16）中把内点p取到边界上P点（三维举例），如图9-3所示。这时，为了避免奇性，需要在P点附近稍加改变，即以P点为中心，以为半径作一小半球包含P点。这样P点仍为内点。此时把边界分为两部分：一部分是鼓起的半球部分，用表示；另一部分即。假如点P取在上（上也一样分析），对于式（9-16）在边上的积分分为两部分，则 （9-17）式（9-17）右端第二工面当时取极限，并把三维基本解代入时，则对光滑表面有当时，而式（9-17）右端第一项仍为上积分，即式（9-16）右端第一项在上积分同样分为两部分，即和。而而 如果P点取在上也可推出同样的结论。所以，把P点取在

8、边界上时，则式（9-16）为 （9-18）上式即泊松方程的边界解公式。它把边界上的函数值与边界上的导数值通过基本解的边界积分之间的关系联系起来了。如果，即为拉普拉斯方程的边界解公式，对于光滑边界，即为 （9-19）或者写成一般式为 （9-20）式中，P，Q皆为边界上的点，C（P）是与P点处的边界几何形状有关的常数，如图9-4所示。图9-4 C（P）与边界几何关系 （9-21）式中的为边界点P处的边界切线之间的夹角。对光滑边界，所以。2边界离散为边界单元并建立代数方程组由上面推得了两个关系式，（i）即式（9-16），得出域内点p的函数值与边界上P点的函数值及导数值之间的关系。换句话说，边界上的函

9、数值和导数值求得之后，即可求出域内任一点的函数值。（ii）即式（9-20），得到边界上任一点的函数值与边界上其它点函数值与导数值之间的关系。下面对式（9-20）进行边界离散化，即划分为边界单元来求边界上的函数值及其导数值。以二维问题为例，离散方案一般有三种情况，如图9-5所示，即图9-5（a）常值单元，（b）线性单元，（c）二次单元。所谓常值单元即在单元中函数值和导数值均为常数，用中个单元具有两个节点，用节点值对对函数作线性插值。二次边界单元，每个单元有三个节点，函数和导数在单元内作二次变化，单元的几何形状适应于曲线边界。下面用常值单元将边界解公式化成节点值的线性代数方程组。如图9-6所示，把

10、边界划分为n个边界单元。其中n1个单元属于；n2个单元属于。单元总数为n = n1 + n2。对于节点Pi，式（9-20）可写为 图9-5 边界元的不同类型 图9-6 常值单元（a）常值单元；（b）线性单元；（c）二次单元 （9-22）由于和q在单元内为常数故可提到积分号外面。Pi为单元上的节点，Q点为单元上的任意点。上的积分与点Pi和单元有关。 （9-23）由于已知（基本解），也可求得，故式（9-23）亦可求得。于是式（9-22）可写成 （9-24）进一步可写成 （9-25）式中 （9-26）把式（9-25）写成矩阵形式，则为式中 （9-27）； （9-28）； （9-29）在式（9-27）

11、中，两边都有未知量和已知量，即在，（边界上函数值给定）；在上，（边界上导数值给定）。给在上有个q值待求，在上有个值待求，把式（9-27）中的已知值和未知值重新整理分别写出，则为已知已知未知 未知 （9-30）把式（9-30）简写为 （9-31）把式（9-31）的未知量移到方程左边，已知量移到方程的右边，则得 （9-32）写成常用的方程组表达式时，即为 （9-33）式中在式（9-33）中X 中包含有个待求的未知函数值和个待求的未知导数值，其未知量仍为n个。F中有个函数值和个导数值作为已知的边界条件给出，n个方程正好求解n个未知数。矩阵A是一个非稀疏的非对称的矩阵，与有限元系数矩阵K是不同的。我们

12、把式（9-16）（这里）改写为离散形式，即 （9-34）一旦整个边界上的值和q值都求得了的话，就可利用上式求得域内任何点p处的值。3系数矩阵H，G中元素的计算首先来计算对角线元素和。是的系数，它是在单元上积分得到的。当时，是边界单元的中点（常值单元），而Q是边界单元上的任意点，如图9-7所示。边界单元的两端点的坐标分别是（x1，y1）和（x2，y2）表示。图9-7 在上的积分路线由式（9-26）及（9-23）第一式知，对光滑边界，有而 （9-35）式中，r为Pi到Q的距离，而n是外表面法线，r和n正交，故，所以于是 是的系数，它也是在上积分的。由式（9-23）知：把拉普拉斯方程基本解式（9-1

13、4）代入，则 （9-36）取元因次局部坐标，在单元中点，在两端点分别±1。令，为单元边界的长度。代入到式（9-36）中，则 （9-37）再来考虑非对角线中和Gij。如图9-8所示。当时，Pi是单元的中点，Q是单元上的任意点，和Gij的计算，可采用数值积分（如高斯四点积分）进行。我们先导出数值积分所采用的公式。由式（9-23）第一式，得图9-8 Pi在上的积分路线由于故 （9-38）图9-9 Pi到的垂直距离把式（9-38）代入式（9-23），得 （9-39）式中，sj为边界单元的长度；rij为Pi到中Q点的距离；hij为Pi点到单元的垂直距离。由图9-9所示，有如果 r和n

14、的夹角为锐角则hij取正号，为钝角则取负号，如图9-10所示。也可由Pi到两端点的矢量r1和r2的叉积来判别，即逆时针时，hij > 0，顺时针时hij < 0，如图9-11所示。 图9-10 hij的正负号 图9-11 由确定hij的正负号如果把r1和r2分别写成 ，则式中，于是hij的正负号可按c值来决定，即c>0，hij > 0；c<0 hij<0。Gij的计算可由式（9-23）第二式得，即 （9-41）由式（9-39）和（9-41）所表达Hij和Gij可用高斯四点数值积分算出，即 （9-42）式中Wk为权系数，如图9-12所表示的积分方案。图9-12

15、 四点数值积分的几何定义在泊松方程中，有一项域积分，如式（9-18）左边第一项。在离散解时，也要把这项并入到式（9-33）右端项F中去。其作法可把求解域划分为单元，然后采用数值积分，即式中，Ae是域单元面积；k从1到7表示采用7点高斯积分方案（见第三章表3-1）；m是域划分的单元数，见图9-13所示。当然，这项积分也可化为边界积分，在弹性全求解时叙述。4域内p点函数值的计算在计算出边界上的值和q值以后 ，利用式（9-34）式可计算域内任一点p的函数值。如果所求为泊松方程，则需加入域积分项，这时式（9-34）为式中，Hij和Gij虽然也由式（9-39）和（9-41）计算，但其中rij要换成域内p

16、点到边界上的Q点的距离。5线性边界单元和高阶单元的计算首先我们分析线性边界单元。线性单元具有两个节点，函数和导数在单元内作线性变化。用单元的相交点取为节点如图9-14所示。 图9-13 域积分 图9-14 线性单元对n个单元，把式（9-20）写成离散形式为 （9-43）由于和q在单元内呈线性变化，不能直接提到积分号外面，故计算要麻烦一些。利用两点插值公式，对和q进行插值，写成下列形式： （9-44）式中，和分别为单元节点j和节点j + 1的函数值及导数值。是无因次单元局部坐标，即如图9-15所示。N1和N2当的函数，分别由下式给出，即图9-15 单元上的坐标下面分别来求式（9-43）的积分。把

17、式（9-44）分别代入式（9-43）中去，左边第二项积分为 （9-45）令； （9-46）于是 （9-47）同理 （9-48）式中 ； （9-49）把式（9-47）和（9-48）代入式（9-43）中去，得 （9-50）上式展开后得 （9-51）对于单连域问题，边界首尾卸接，这时； （9-52）将式（9-51）整理后得 （9-53）令 ； （9-54）式中 ； 于是，对于节点Pi，式（9-43）可写成 （9-55）与常值单元一样，进一步可写成移项后，得AX = F （9-56）下面导出微元的计算，式（9-46）和（9-49）N1和N2是的函数而积分微元是对进行的，故必须把变换到局部坐标系中去，如

18、图9-16所示。与总体坐标（x，y）的关系为图9-16 曲线边界的几何定义由于所以 （9-57）称为雅可比行列式。如果对总体坐标x，y也进行线性插值，即 （9-58）式中，（xj，yj）和（xj+1yj+1）分别为边界单元两端点的坐标值。由式（9-58）可求得代入式（9-57），得 （9-59）式中sj为单元的长度。于是，对于线性单元，式（9-46）和（9-49）分别由下式表达， 即 （9-60）式（9-60）中的积分同样可用四点高斯积分计算。图9-17 二次单元二次单元是在边界单元内的函数和其导数均采用二次插值，如图9-17所示，其式为 （9-61）式中，为的二次函数，即 （9-62）与线性

19、单元类似的作法，求下面的积分： （9-63）式中 （9-64）设，分别为单元1，2，3节点的总体坐标。与线性单元一样，可写成 （9-65）于是有 （9-66）将式（9-66）代入式（9-57），得 （9-67）其它推导与线性单元相同，但在每个单元的首尾节点（如节点1和2），同样有两项系数叠加以组成总体方程组的系数。三、弹性体受力分析的边界元法1基本关系和预备知识为了简化起见，在弹性力学的表达式中用张量表示。所以，首先把用到的张量符号简介如下。指标记法：在直角坐标系中，通常用x，y，z表示坐标。 如果分别用x1，x2，x3来表示上述三个坐标，记为xi（i =1，2，3）就称指标记法。这里i称为下

20、标。ai（i =1，2，N）代表a1，a2，aN N个量，则代表N×M个量，则代表N×M×L个量，依此类推。求和约定及哑标：某指标在某一项中重复出现，且仅重复一次，则该项代表一个和式，按重复指标的取值范围求和，如表示：表示求和重复的指标为哑标，用什么字母均可。不求和的指标称为自由标。如上例i就是自由标。克罗内克符号： 定义为=1（i = j时）；=0（i j时）当i, j =1，2，3时，表示可知； ； ； （当i, j, k为1，2，3循环序列）（当i, j, k为1，2，3逆循环序列）（当i, j, k有相同指标）置换符号eijk定义如下：eijk共代表27个

21、量，有21个量为零。显然：导数记号：求导符号用“，”表示。中用指标记法表示弹性体的基本关系时，其平衡方程记为 （在内） （9-68）力的边界条件为 （在上） （9-69）位移边界条件为 （在上） （9-70）几何方程为 （9-71）各向几性弹性体的物理方程为 （9-72）式中 ； （9-73）称为拉密常数；G称剪切模量。或者更简明地写为式中 用位移表示的平衡方程为 （9-74）如果考虑热影响时，则用位移表示的热弹性方程为 （9-75）考虑热影响时的应力表达式为 （9-76）式中 式中，a是热膨胀系数，T是温差。以上指标记法，对于三维问题，3；对于二维问题，i，j，k=1，2。对于平面应力问题需

22、将弹性常数作相应改变，即使E，a相应变为，。2弹性问题边界积分方程表达式如前所描述的势问题一样，我们把平衡方程和两类边界条件写成加权余量表达式为 （9-77）式中，和为权函数，分别为加权场的位移和面力，也即基本解。同样有 （9-78）现在把式（9-71）和（9-72）应用到近似解和中权场中去。对式（9-77）左边第一项进行分部积分，便得 （9-79）如果考虑温度影响时，则需把式（9-75）代入到式（9-79）中左边第一项中去，即 （9-80）有关温度影响，可并入体积力的域积分，在后面推导体积力时详细论述，在此省略。对式（9-79）再次进行分部积分，且用互换定理 （9-81）时，便得到 （9-8

23、2）方程左边体积力的域积分不含未知函数，有关它的计算在以后介绍。现在的问题是如何把左边第一项转变为边界积分。弹性问题的基本解是满足平衡方程 （9-83）的解。这里 函数是表示在i点在l方向施加的单位载荷。如果对上述作如下运算，即于是得 把上式代入式（9-82），得 （9-84）式中，表示域内任意点i（或p）沿l方向的位移；为基本解，即在无限域中式（9-83）的解。其意义分别为在i点沿l方向有单位力作用时所产生的位移和表面力。如果考虑i点的三个方向时，则式（9-84）为 （9-85）或者把i换p表示时，也可将上式写成 （9-86）式中，和分别在域内p点沿l方向作用单位力时，而在边界上的Q点在k方

24、向产生的位移分量和面力分量。l，k分别代表p点（或i点）和Q点的有关坐标方向。对三维问题，l，k = 1，2，3。l，k的符号意义如图9-18所示。对二维问题，l，k = 1，2。式（9-86）是弹性体的边界积分方程。该式表明区域内任意p点的位移分量ul（p）和边界上Q点的位移分量u（Q）和表面力分量pk（Q）之间的关系。如果边界上所有未知的位移分量和表面力分量已求得，那么域内任意点的位移分量可以通过该式求得。图9-18 在p点在l方向（l = 1）作用单位力而在Q点所产生的3二维弹性体的基本解和边界元方程的表达式二维弹性平面问题的基本解，可由无限域的式（9-83）导出。对于平面应变问题，位移

25、和面力的基本解分别为 （9-87）对于平面应力问题，则分别为 （9-88）式中，r是单位载荷作用点p到所考虑的点Q之间的距离。nj是该点的表面法线的方向余弦。n为外表面法线。当l = k时，时，。下面来推导边界点的积分方程。式（9-86）对区域内任意点都适用。当p点取在边界上P时，基本解会产生奇性。与处理势问题一样，以奇点P为圆心以为半径作小半圆，如图9-19所示。于是边界又分为两部分；一部分是；另一部分是。把基本解代入式（9-86），当取极限时，式（9-86）右边第一项在上积分为 零，而在上积分仍为上的积分。式（9-86）左边第二项在上积分，当 它与边界形状a角有关。积分后代入到式（9-86

26、）中去，得 （9-89）图9-19 p点取到边界P时的处理式中，P即所取的边界点；而为边界上的位称移；Clk是与边界几何形状有关的常数，表达式为 （9-90）式中，对于光滑边界，所以 （9-91）式（9-89）是边界上的位移分量uk（P），uk（Q）和表面力分量pk（Q）之间的边界积分关系式。由该式可以求出边界上的全部未知位移分量和表面力分量。与势问题一样，我们把式（9-89）的边界积分进行离散化，可离散成常值单元、线性单元及二次和高次单元。该式写成矩阵表达式为 （9-92）是的2×2矩阵，P*是的2×2矩阵，即； （9-93）； ； （9-94）下面我们把式（9-92）离

27、散成n个单元，并采用常值单元，则该式为 （9-95）式中，Ui为边界单元的位移矢量，Uj，Pj分别为边界单元的位移矢量和面力矢量。由于采用的是常值单元，和势问题一样节点取在单元中点，位移变量和力变量在单元内是常值，故可提到积分号外面去。n是边界单元数。右端第二项为体积力的域积分。令 （9-96）式中，m为域积分单元数，l为域积分的积分点数，Wk为权系数，Ae为单元e的面积。这是把体积力的域积分用加权数值积分得出的。当然，也可转变为边界积分计算，在后面专门介绍。于是式（9-95）可写成 （9-97）进一步可写成 （9-98）式中，当时，；当i = j时，对光滑表面，方程（9-97）或（9-98）

28、对节点i给出一组方程式（三维是三个，二维是两个），Uj和Pj都是在节点j上的未知数。和Gij是节点i与物体表面上全部节点相关的系数。这里与势问题的区别在于每个节点位移和面力有两个分量（三维问题是三个），故和Gij是2×2子矩阵。对所考察的每个节点都写出式（9-97），于是有 （9-99）式（9-99）写成简化形式为 （9-100）求解方程组（9-100）以前，必须加入边界条件，此处的边界条件为在上 （）在上 （）式（9-100）共有2n个位移值和2n个面力值，但有n1个位移值和n2个面力值是作为边界条件给定的，因而式（9-100）中尚有2n个未知量，即有n1个面力值和n2个位移值待求

29、，把已知位移和已知面力代入式（9-100）中去，使方程组重新组合，即将所有未知项表示矢量X，写在等式的左边，其结果为 （9-101）体力矢量项也包括在X中。4系数矩阵中元素Hij和Gij的计算首先计算对角线元素Hij和Gii。当j = i时，式（9-96）中有关边界积分是在的中点，而Q点取在单元中的任意点（在同一个单元中），如图9-20所示。图9-20 Pi和Q在中积分由式（9-93）及（9-96）可知 （9-102）把式（9-88）代入上式，考虑到，同时由图9-20可见，在中，。在中，。代入计算后，便得到。于是得而ci也可由研究刚体运动求得。若设在任一方向上刚体有一个单位的位移，则式（9-1

30、00）变为式中，Il是在l方向上刚体位移的单位矢量，于是H的 对角线项可简写为这意味着并不需要用显式来确定系数，对光滑表面。Gii的计算是由式（9-93）和（9-96）中第二式进行的，即 （9-103）把式（9-88）中的代入上式积分，注意，当r = 0时为奇点，故积分限可取为到，使时取极限，计算后得 （9-104）下面介绍非对角线元素Hij和Gij的计算。当时，Pi是取在边界单元上的中点，而是在其它边界单元上的任意点，如图9-8所示。计算式仍由式（9-102）和（9-103）计算，只是把下标ii换成ij，也是2×2矩阵。把式（9-98）中的和分别代入和Gij中，计算中令，（见图9-

31、8），于是具体表达式分别为 （9-105）和 （9-106）假如采用4点高斯积分计算上面积分，只需将微元代入上式，如其它，类似计算。sj为单元边界长度；Wk为权系数；h为Pi点到单元垂直距离；r为Pi到单元积分点的距离；n1，n2分别为的方向余弦。5域内点的位移和内点、边界点应力的计算一旦边界上的位移和面力算出后，就可以用式（9-86）求得域内点的位移。把该式写成离散形式时，则为 （9-107）式中，Ui为内点Pi的位移矢量。Hij和Gij分别由式（9-105）和（9-106）计算，但式中r为内点Pi到 边界点之间的距离。下面引出内点应力和边界点应力的计算公式。由式（9-72）可知，应力表达式

32、为 （9-108）将式（9-106）代入几何方程（9-71），再代入物理方程（9-72），得内点应力表达式为 （9-109）上式是全部在内点进行的。令 （9-110）于是应力表达式为 （9-111）3阶张量Dkij和Skij的推导简要说明如下。首先将在p点对xj微分： （9-112）求得 ； 代入式（9-112），于是得 （9-113）式中，如果将下标j换为l，且，则 （9-114）将式（9-113）的下标l换为i，便得 （9-115） 如将式（9-113）的下标l和j分别换为j和i，就可以得到 （9-116）将式（9-114），（9-115），（9-116）代入式（9-110）整理合并后，便

33、得 （9-117）现在推导Skij表达式。在p点将对xj微分得 （9-118）将上面关系以及前面导出，等代入式（9-118），得 （9-119）将上式中下标j换为l，便得到；将下标l换成i，便得到；将下标l和j分别换成j和i，便得到。将其结果均代入式（9-10）第二式中去，得 （9-120）当采用常值单元并把式（9-111）写成离散后的矩阵形工，得任意内点p的应力为 （9-121）式中 ； ；中的系数为，；中的系数为，；中的系数为，。它们分别表示如下： （9-122）和 （9-123）式中r为内点p到边界点Q之间的距离。上述公式在编程序计算时，可用4点高斯积分解决，例如：计算和可把式（9-12

34、2）和（9-123）的第一式分别写成：和其它类推，只是把积分写成和式，把而已。在数值积分式中Wk为权系数。边界上的应力在工程结构中往往是感兴趣的。边界上的全部位移和面力求出后，就可以计算边界上的应力。由于边界上的位移uk（Q）沿边界s上的变化是可以求得的（近似求出可）。再加上物理方程（平面问题为三个） 和力的边界条件（二个），七个方程正好求出七个未知数，即三个应力和四个位移在边界上的导数。具体写为 （9-124）式中，l = 1，2。sl是边界的参数，在上是已知的（差分计算或对uk直接微分）。把上式写成矩阵形式为 （9-125）式中，；是垂直s的单位矢量；s是弧长。四、体积力的边界积分计算 体

35、积力的计算，在机械工程中是经常遇到的。在热动力机械中，体积力包含有重力、离心力和温度载荷。含有体积力的弹性体的边界积分方程是式（9-86）或（9-89）。由该二式看出，有一项是域积分。对这项的处理，如果采用在域内划分单元，然后采用数值积分如式（96-96）那样，则是不胜其烦的，且所需的数据准备也较多。如果体积力为常量或无热源的稳态热载荷，就可以把体积力的域积分转化为边界积分。这就充分利用了边界单元法的优点。现在我们先针对二维问题来讨论这种算法。域积分的积分项在积分方程中为 （9-126）式中，是整个弹性体的体积域。为了把它转变成边界积分，我们在这里定义一个张量Gki（P，Q）p是力点，Q是场点

36、。这个张量叫伽辽金张量。二维问题的表达式为 （9-127）由它组成的基本解为 （9-128）把式（9-127）代入到式（9-128）中去，得 （9-129）可以看出式（9-129）与（9-88）只差一个常数。常数只表示刚体位移，对计算结果是无影响的。把式（9-128）代入到式（9-126）中去，得 （9-130）这个公式似乎复杂些，但易于转变为边界积分。下面我们针对几种常见的体积力把式（9-130）转化为边界积分公式。1重力载荷在常重力场中，由于为常数，为质量密度，为加速度。根据高斯定理，有 （9-131）这就是重力为体积力的体力项转变为边界积分的表达式。把式（9-127）微分后代入式（9-1

37、31）中去，便得 （9-132）于是，内点边界积分方程可以写为 （9-133）式中 （9-134）在重力场作用下的应力表达式为 （9-135）式（9-135）中的是这样求得的，首先把式（9-134）代入到应力位移关系（即式（9-72）中去，得到 （9-136）再把式（9-134）求导代入式（9-136）得到 （9-137）2离心力载荷在二维平面情况下，旋转体单位体积所产生的离心力为 （9-138）令 =常数把式（9-138）代入到式（9-130）中去，得 （9-139）现在对式（9-159）作如下变换：二式相加得到等代入式（9-139），便得到 （9-140）用高斯定理把域积分转变为边界积分，

38、即 （9-141）由于是对称的，方程（9-141）稍加改变，便得 （9-142）如重力载荷情形一样，民可把式（9-142）简写为 （9-143）式中 （9-144）应力表达如同式（9-135），对于离心力所产生的则为 （9-145）3热载荷对二维稳态热载荷，体积域积分是由温差所产生的。根据热弹性方程的马克斯威尔-贝蒂互换定理，得内点边界积分方程为 （9-146）于是 （9-147）式中，由式（9-128）得 （9-148）由改变整数虚标，于是式（9-148）变为 （9-149）代入式（9-147），得 （9-150）对稳态热载荷，有 （9-151）于是，我们可以把式（9-150）写成 （9-1

39、52）由高斯定理，可以把上式转变为边界积分，即 （9-153）这就是由于温度载荷所产生的边界积分式。该式也可简写为 （9-154）式中 （9-155）把式（9-127）代入式（9-155），分别得 （9-156）于是，对于热弹性方程位移边界积分表达式为 （9-157）把式（9-157）微分后代入应力位移关系中去，便得 （9-158）式中 （9-159）五、三维弹性体的边界元公式1三维弹性体的基本解和三阶张量弹性力学三维问题的基本方程为式（9-68）（9-76）诸式，只是其中下标的取值范围分别为。积分方程表达式的推导与二维问题是类似的，在此不作重复。域内任意点p的位移ul（p）与边界积分关系以及

40、边界点位移cul（P）与边界积分关系仍为式（9-86）和（9-89），只是其中l，k =1，2，3。式中的基本解，则分别为 （9-160）（9-161）式中，见图9-21所示。图9-21 几何意义对于光滑表面，系数c仍为，非光滑表面可由二维推广。其它如边界方程离散化，方程组的建立，与二维问题大体是类似的。只是现在自由度增加了，每个边界节点有三个位移分量和三个面力分量，n个节点有3×2n个方程式。积分边界是三维空间域表面，即面积分。可采用有限元法有关插值公式和数值积分方案进行。在式（9-111）应力表达式中的三阶张量Dkij和Skij分别为 （9-162）式中，导数取在边界上，。2三维

41、体积力的边界积分式三维体积力边界积分公式的导出，与二维是类似的。相应于三维基本解的伽辽金张量为 （9-163）将上式微分两次代入式（9-128），得 （9-164）重力载荷的三维体积力边界积分式：把式（9-163）微分两次代入式（9-130）中，得 （9-165）式中 （9-166）用与二维问题类似的方法得应力表达式中的为 （9-167）离心力载荷三维体积力边界积分式：与二维问题一样，离心力体积力边界积分也可写成 （9-168）式中 （9-169）应力表达式中的为 （9-170）三维稳态热载荷的边界积分式：在稳态热载荷中，体积力的边界积分项也可写为式（9-154）的表达式，即式中 （9-171

42、）应力项中的和分别为 （9-172）六、轴对称体的边界元公式燃气涡轮发动机等机械含有大量的轴对称体零、部件。这些零、部件用边界元法分析其结构完整性，会使其输入数据和计算时间大为减少。1轴对称体基本方程及伽辽金矢量用圆柱坐标表示的轴对体的纳维叶方程为 （9-173）式中其矢量形式为 （9-174）为了求解方程组（9-173），定义伽辽金矢量为 （9-175）把式（9-175）代入式（9-173）中，得双调和方程为 （9-176）式（9-176）写成标量时，则为 （9-177）式中，Gr和Gz是伽辽金矢量的径向和轴向分量；Fr和Fz是体积力分量。这里的体积力项可通过应用Dirac 函数来表示，即

43、（9-178）式中，R和Z是载荷点的坐标；r和z是场点的坐标。根据Dirac 函数的性质，有 （当时）为了得到边界解的公式，需要找到轴对称体纳维叶方程的基本解。由前述，对一般三维问题，其基本解是均匀无限体的点载荷的开尔文问题的解。对于轴对称问题，基本解的物理意义是轴对称环载荷的开尔文解。今考虑一个具有任意横截面形状的轴对称体，在体内p点为载荷点，其坐标是。边界点为Q点，其坐标为），如图9-22所示。小写p表示固定的载荷点，大写Q表示可改变位置的边界点。其基本解是描述Q点的状态与p点的环载荷之间的关系。图9-22 轴对称体有两种方法可以求得轴对称体的基本解。第一种方法是围绕对称轴沿环载荷路径积分

44、在维的基本解；第二种方法是直接解出轴对称形式。这里介绍的是第二种方法。用积分变换方法，可以证明，下列解 （9-179）是满足无体力情况下方程（9-177）的。式中称为第二类 阶v次的勒让德函数。根据递推关系可以用零阶的勒让德函数来表示-1阶的勒让德函数，即 （9-180）如果把式（9-180）代入式（9-179）中去，将其结果代入式（9-175）便得到用零阶勒让德函数表示的轴对称体位移核函数。在具体计算时，可用第一类和第二类椭圆积分来表示勒让德函数及其导数。2轴对称弹性问题的基本解（位移核和力核）为了方便起见，用完全椭圆积分来表示伽辽金矢量而不用勒让德函数时，便得到 （9-181）式中，和分别是第一类和第二类完全椭圆积分并定义为 （9-182）式中把伽辽金矢量分量代入式（9-175）中，得位移矢量为 （9-183）式中，er和ez分别是径向和轴向的单位矢量。在边界Q点，径向和轴向的位移分量为 （9-184）这里pr和pz分别是内