边界元法

边界元法边界元法是继有限元法之后的一种独具特色的新型数值方法，它是将描述数学物理问题，包括力学问题的偏微分方程边值问题，化为边界积分方程，并吸收了有限元法的离散化技术而发展起来的。边界元法中包含有有限元法的思想，它将有限元法的按求解域划分单元离散的概念移植到边界积分方程的方法中，但边界元法不是有限元法的改进或发展，边界元法与有限元法存在着质的差异。有限元法是在整个求解域上进行离散，边界元法只在求解域的边界上进行离散；有限元法是全域数值方法，而边界元法在域内采用了物理问题或弹性力学的基本解和一些积分运算，数值计算只在边界上进行，它属于半解析半数值方法。由于边界元法具有的一些有限元法所没有的优点，人们对其的研究兴趣和信心不断增强，研究的范围也不断扩大。边界元法已经被用来求解各种力学和非力学问题、线性和非线性问题等。§1边界元法概述Rizzo于1967年提出对弹性力学边界积分方程按边界离散的思想，是弹性力学边界元法研究的开端。尽管有限元法所取得的成就与日俱增，但有限元法还不是十全十美的。改进有限元法的努力一直在进行着，但有限元法的某些不足是无法克服的。例如有限元法需全域离散，导致问题的自由度和原始信息量大；对无限域只能人为地取成有限域；有限元法的离散技术本身也存在缺陷，它把本来是连续的介质用仅在节点处连接的有限单元的集合来模拟，这样不仅带进了离散的误差，而且在单元之间连续的要求较高时，有限单元的构造也往往很困难；对同一问题采用不同的程序计算时可能会得出不同的结果，对有限元法的精度和可靠性也常常会提出疑问。相对于有限元法这些不足之处，边界元法仅在边界上离散，使数值计算的维数降低一维，从而减少了问题的自由度和原始信息量。边界元法采用无限域的基本解，用边界元法求解无限域问题可称是天衣无缝。边界元法的离散误差只产生在边界上，边界元法中部分采用数值法，部分采用解析式，它的准确度和可靠性已公认是高于有限元法。此外，边界元法在弹性力学中扎的根比有限元法还要深一些。边界元法的研究在经过对弹性力学和塑性力学问题的初步尝试后，没有能够得到满意的计算结果。边界积分方程和从有限元法借鉴来的边界离散技术两者不是一试就灵的。较为完整的、可以实际应用的边界元法是20世纪70年代才建立起来的。先后解决了边界积分方程中含有的奇异积分的精度问题，对于不带奇异性的积分，只要用普通的高斯数值积分法就可以获得足够的精确度，从而使边界元法固有的精度得以显露出来。20世纪80年代可称是边界元法蓬勃发展和丰收年代。边界元法的发展大致可归结为以下三个方面：（一）数学方面。包括边界元法的数学分析理论和数值积分方法的研究。边界元法的发展固然是由于电子计算机的迅速发展和广泛应用使其实施成为可能，也与近代数学理论的发展密切相关。边界元法数学方面的研究，不仅克服了由于积分奇异性造成的困难，同时又对收敛性、误差分析以及各种不同的边界元法形式的统一进行了数学分析，为边界元法的可行性和可靠性提供了理论基础。边界元法数学理论方面的主要工作有：边界元法包括数值积分的渐近误差分析，有限元和边界元耦合方法的误差分析，边界元法解的稳定性，弹性和流体力学问题的边界元法的数学理论，断裂力学边界元法的误差分析等。近年来，一些学者将有限元法理论中的区域分解方法引入边界元法，讨论了其数学理论。总的说来，边界元法数学理论的研究还落后于方法和应用的研究，与有限元法数学理论的研究尚有一定的差距，有待进一步和发展。（二）方法应用方面。包括边界元方法的完善和应用范围的拓宽。20世纪70年代以前，边界元法的研究只限于解决以下几个方面的问题：势问题、弹性静力学、瞬态弹性动力学、波的传播、断裂力学、流体力学、板弯曲问题等，而且对一些问题的研究只是初步尝试。现在，边界元法的发展已涉及工程科学的很多领域，边界元法的应用也已经规范化。对非线性问题，其方法亦趋于成熟。在工程和工业技术领域，边界元法的应用已涉及到水工、土建、公路、桥梁、机械、电力、地震、汽车、航空航天等诸多方面。（三）应用软件方面。边界元法作为一种数值计算方法，其应用要通过计算程序来实现。这种计算程序作为应用软件，是随着边界元方法的发展而发展的。20世纪70年代以后，随着边界元法国际会议在世界各地逐年举行，陆续有边界元法应用软件的新成果问世。1982年第四届边界元法国际会议上，英国南安普敦大学的Danson介绍了他们研制的边界元分析程序包BEASY，它是国际上第一个边界元的大型软件。1985年以来，边界元技术国际会议在世界各地举行，它着重于边界元计算技术的研究和应用，包括工程应用、计算技术和工业应用等，这对边界元应用软件的发展起到促进作用。现在，边界元法应用软件已由原来的解决单一问题的计算程序向具有前后处理功能、解决多种问题的边界元程序包发展。已经形成的程序包有BEASY(英国)、CA.ST.OR(法国)、BETSY(德国)、SURFES(日本)、ESBEA(美国)等。但是，目前边界元程序包无论在质量上还是数量上均与有限元法程序包的发展有一定的差距。要想使边界元法像有限元法那样得到广泛应用，还必须发展各种各样的高质量的通用程序包。我国关于边界元法的全面研究始于1978年。几十年来，随着国际边界元法的发展，我国学者也取得了很大成绩。在数学方面，我国开展的主要工作有：自然边界元法的数学理论；椭圆型方程边界元法的数学理论；边界元法小区域分解方法；收敛性和误差分析等。在边界元的方法与应用方面，我国学者在求解各种问题的边界元法的研究方面做了很多工作。主要有：弹性力学、弹性动力学、断裂力学、动态断裂力学、板壳问题、势问题、流体力学、非弹性力学、弹塑性有限变形、波的传播、复合材料力学、复合材料界面断裂力学、岩土力学有限元与边界元耦合方法和反问题等。在应用软件方面，伴随着边界元法的研究，对以上求解各种问题的边界元法，我国学者都发展了应用软件，有些已经应用于工程实际问题，并取得良好的效果。但是这些软件多是针对某一类特定问题，没有具备完善的前后处理程序，因而应用上有些还不如有限元法通用程序包方便。边界元法有直接法和间接法两种。直接法通常是用格林恒等式或加权残值理论来表述，采用的变量物理意义明确，并能通过边界积分方程数值离散后直接求解，是边界元法的主要方法。间接法则是利用位势理论来推导公式，使用的变量物理意义不太清楚。当然，间接法仍有它的可取优点。本章将介绍这两种方法的主要思想和实现过程。§2直接边界元法本节将给出势论问题的直接边界元方程的建立和推导过程，对于其他的工程问题，基本概念都是相同的。考虑如图2.1所示的位势问题，区域为Ω，边界为21ΓΓ=Γ∪。图2.1位势问题势函数φ所满足的基本控制方程（Laplace方程）以及边界条件表示如下02=∇φ在Ω中（2.1）φφ=在1Γ上（2.2a）qn=∂∂φ在2Γ上（2.2b）在本书的第三章介绍了加权残值法的概念和应用。这里，我们利用加权残值法来推导边界元方程。首先引入权函数*φ，其具有连续的一阶导数，则加权残值的提法表示为Γ∂∂−−Γ−∂∂=∇∫∫∫ΓΓΩd)(d)(d)(***212nqnsφφφφφφφ（2.3a）上式中昀后一项是用边界1Γ上的条件加权得到的。对上式中的关于区域Ω内的积分项，进行两次分部积分运算后，得到Ω)(1φφ=Γ)(2qn=∂∂ΓφxzonsyΓ−Γ∂∂+Γ∂∂−Γ∂∂=∇∫∫∫∫∫ΓΓΓΓΩddddd)(*****22211φφφφφφφφφqnnns（2.3b）由上式可以看出，只要找到满足一定条件的权函数*φ，上式中的区域积分项可望获得进一步的简化。这里，我们寻找在区域Ω中满足控制方程0*2=+∇iδφ（2.4）的解，这个解称为基本解。ir和r分别表示位势问题中的源点i处的矢径和场点矢径，)(iirr−=δδ表示狄拉克函数。对于三维问题，算子rrr∂∂+∂∂=∇2222，不难得到（2.4）的解为rπφ41*=（2.5）对于二维问题，算子rrr∂∂+∂∂=∇222，其对应的解为)1ln(21*rπφ=（2.6）将（2.4）代入（2.3b）式中的左端，我们可以得到isφφφ−=∇∫Ωd)(*2（2.7）从而（2.3b）式可写成Γ+Γ∂∂=Γ∂∂+Γ∂∂+∫∫∫∫ΓΓΓΓdddd****2112φφφφφφφφqnnni（2.8a）或者∫∫ΓΓΓ∂∂−Γ∂∂=**ddnniφφφφφ（2.8b）其中Γ表示区域整个边界，并在1Γ上有φφ=,在2Γ上有qn=∂∂φ。上式表明区域中任一点的位势值可用边界上的φ和n∂∂φ来表示。在上面的方程（2.8）中，如果能将i点移到边界上去，那么方程中的全部量都是边界上的量，这也正是建立边界积分方程的基本思想。如图2.2所示，将点i移到边界处，但点i附近的边界要改变一下，变为以点i为球心的鼓起的部分球面（以三维情形为例）。球的半径r很小ε=r。此时，我们将边界分为两部分，一部分是鼓起的球面εΓ，另一部分是剩下边界Γ′。niΓiΓ′iθεεΓ图2.2边界上的处理对于新的边界Γ′Γ∪ε来讲，点i依然为内点，（2.8）式仍可以用，即∫∫∫∫Γ′Γ→→Γ→Γ′→Γ∂∂−Γ∂∂−Γ∂∂+Γ∂∂=*0*0*0*0dlimdlimdlimdlimεεφφφφφφφφφεεεεnnnni（2.9）这里需要解决上面式子中的边界εΓ上的积分0d4limd41limdlim00*0=∂∂=Γ∂∂=Γ∂∂∫∫∫Γ→Γ→Γ→εεεθπεφπεφφφεεεnnn（2.10a）iiuπθπεφπεεφεεεε4d41limd)41(lim200−=Γ−=Γ∂∂∫∫Γ→Γ→（2.10b）上式中利用了θεdd2=Γ，iθ表示鼓起部分球面对点i所张的立体角。当0→ε时，部分球面收缩于点i时，边界Γ′趋向于原来的Γ。对于二维问题，可以类似地进行处理。昀后，（2.9）可以表示为∫∫ΓΓΓ∂∂−Γ∂∂=**ddnnciiφφφφφ（2.11）并且⎪⎪⎩⎪⎪⎨⎧−−=):3(41):2(21为立体角为平面角iiiiiDDcθπθθπθ（2.12）对于光滑边界，系数2/1=ic将变得更为简单。显然（2.11）给出了边界Γ上φ所满足的方程，方程中全部的量均为边界上的物理量。§3间接边界元法在边界解法中，昀简单的一种解法即为间接法。这种方法简单而且适用，对于许多实际问题能够给出较好的结果。但直接法与之比较，则更为一般，应用范围更加广泛。间接法可以看作是直接法的特殊情形，解释为加权残值法的特殊情形。这一方法的具体实现，首先从满足区域中控制方程但尚有一些待定系数的解入手，然后由满足控制方程的近似解在边界上一定数目的点上通过与边界条件相匹配来决定这些待定系数。对于势论问题，一般由边界条件φφ=在1Γ上（3.1a）qn=∂∂φ在2Γ上（3.1b）来确定待定系数。我们记区域Ω之外的区域为−Ω，在区域−Ω内没有源点。设−φ为−Ω内满足Laplace方程的解，即02=∇−φ。间接法可以由直接法边界方程（2.8b）导出。由于−Ω相对于区域Ω为外部区域，所以当点i在区域Ω内时，对于区域−Ω有0=c，从而经过类似与直接法的推导，可以得到0dd**=Γ∂∂+Γ∂∂∫∫Γ−Γ−nnφφφφ（3.2）注意到上式中的第一项符号与（2.8b）相反，这是由于n∂∂−φ是对区域Ω的内法线求导，而n∂∂*φ是对区域Ω的外法线求导。下面我们分两种情况来讨论：（1）内、外区域的位势函数值在边界上连续。也就是说，区域Ω关于位势问题的边值φ，与区域−Ω关于位势问题的边值−φ相等，即−=φφ。将（2.8b）与（3.2）相加，得到对于区域Ω内部的点i成立的关系式∫Γ−Γ∂∂+∂∂=*d)(φφφφnni（3.3a）如果令nn∂∂+∂∂=−φφω，则上式另表示为∫ΓΓ=*dωφφi（3.3b）上式即为间接边界元法提法的表达形式，内点位势值通过边界量来表达。其中的ω是边界Γ上*φ的待定源密度分布，在用上式确定iφ时它是必不可少的，并且必须先求出来。ω的物理意义是内部和外部问题的解在边界上的导数之差。值得注意的是，无论是对Ω内的*φ还是对于−Ω内的*φ，它们都是由Ω内的点i处点源产