有限元与边界元(下)

第二部分边界单元法第一章边界元法数学基础1.1狄拉克δ函数在物理现象中我们用到的质点，点电荷等物理概念。如质点，为质量集中的点，其体积趋于零，故它的密度（质量/体积）趋与无穷大，但它的密度的体积分（总质量）却是个常数，所以可以任意选择该常数为1。于是得一个单位质点；在电学中的点电荷也有相似的特征，即其体积趋于零，故它的电荷密度g（电量/体积）趋于无穷大，但它的电荷密度的体积分也是一个常数，也取为1。于是得到一个单位点电荷。为了描述这一类抽象概念。在数学上引入了一个δ函数，定义δ函数如下：设)(zyx、、及)(00000zyx、、是区域内的任意两点，0是一个固定点。如果000,0,)(　　　　　　（1.1.1）000,0,1)(MMd　　　　　　(1.1.2)则称)(0为函数。可以把)()()()(0000zzyyxx视为质点或点电荷的位置坐标，那么函数)(0反映了上述物理现象。对于二维、三维函数可以写成下列形式：二维)()()(000yyxx三维)()()()(0000zzyyxx由上述定义式（1.1.1）和(1.1.2)可推出今后我们要用到的几个重要性质。1、当p在域的边界上，且P处的边界光滑时，则21)(dp(1.1.3)证明：对于二维情况，域为一平面，如图1.1.1设P位于域的边界上，我们以P为圆心，在内作一半径无限小的半圆（因P外边界光滑，故内，是半圆）。由（1.1.1）式，)(p在以外都为0，故Pdpdp)()(由(1.1.2)式，因为为半圆，故21)(dp2、设函数u在P点处连续，当P在域内时，则)()(pudpu(1.1.4)式中u(p)为P处的函数值。证明：以P为圆心，在内作一半径无限小的区域包围P点，则由(1.1.1)式有图1.1.1P点位于边界上dpudpu)()(因为为无限小，且u为连续函数，故可用P点的函数值u(p)代替域中的函数值，并移至积分号外。根据函数的定义式(1.1.2),有)()()()(pudppudpu3、设函数u在P点处连续，当P在的边上，且P处边界光滑，则)(21)(pudpu(1.1.5)将性质1、性质2结合起来，很容易证得性质3。1.2格林公式如果u和在区域()上连续且一阶连续可微，在区域内二阶连续可微，则下述格林第二公式成立dnunuduu)()(22(1.2.1)式中为区域的边界（对于三维区域，为边界面；二维区域，为边界线）；n为边界的外法向。证明：由哈密顿算子的运算规则有uuu2)(由区域积分与边界积分的关系有dnudu)(而uuu2)(dnudu)(两式相减，得dnunuduu)()(221.3基本解基本解定义：微分算子L对函数u进行某种微分运算，构成或一个微分方程：L(u)=0(1.3.1)若某函数经过L的微分运算，得到一个()p函数，即()()Lp(1.3.2)则称为L(u)=0微分方程的基本解。微分方程的基本解不是唯一的。例如，如果iu是(1.3.1)式的任意解，而*u是式(1.3.1)的基本解，则)())((0)()()(**ppuLuLuuLii也就是说)(*uui也是(1.3.1)的基本解。下面给出常用的几种微分方程的基本解，它们在处理许多地球物理勘探问题时，是非常有用的。1、三维拉普拉斯方程的基本解三维拉普拉斯方程02u(1.3.3)其中是三维哈密顿算子，它的基本解为r41(1.3.4)其中r是三维区域中某点P至中任意点的距离（图1.3.1）。证明：当0r时，将直接带入球坐标中的2表达式：)(1222rrrrrp)41(1222rrrr=0得2104r当r=0时，r1奇异，上式不成立。我们作一个半径无限小的球面包围P点，2的积分为：ddrd)41(41)41(22根据区域积分与边界积分的关系有图1.3.1三维区域示意图dnrrdrddr)1(41)1(41222114144drrr式中1nr。（由于在域中nr）。可见)()41(2pr(1.3.5)即r41是三维拉普拉斯方程02u的基本解。2、二维拉普拉斯方程的基本解是r1ln21(1.3.6)证明当0r时，将直接带入极坐标中的2表达式：)(12rrrr得0)21(1)1ln21(22rrrrrr当r=0时，r1ln奇异。上式不成立，我们作一个半径无限小的圆包围p点，2的积分为drd)1ln21(22dr)1ln(21根据区域积分和边界积分的关系有：dnrrdrddr)1(ln21)1ln(21212122rrdrr可见)()1ln21(2pr即r1ln21是二维拉普拉斯方程02u的基本解。3、二维赫姆霍兹方程的基本解二维赫姆霍兹方程为022uku(1.3.7)它的基本解分别为)(410krN(1.3.8)01()2Kkr(1.3.9)其中2k为正实数；为(1.3.7)式取“+”号十的基本解；为(1.3.7)式取：“-”号时的基本解；0N是第二来零阶贝塞尔函数；0K是第二类零阶修正贝塞尔函数。4、三维赫姆霍兹方程022uku的基本解是4ikrer(1.3.10)1.4高斯积分法考虑区间[-1,1]上的求积公式:111)()(nkkKxfAdxxf(1.4.1)由勒让德多项式的理论可知,n阶勒让德多项式nnnnnxdxdnxP)1(21)(2！(1.4.2)是[-1,1]上的正交多项式,它在[-1,1]内有n个根nxxx,...,,21,取它们为(1.4.1)的求积节点.则(1.4.1)称为高斯——勒让德求积公式,并且可以推出高斯——勒让德求积系数22)]()[1(2knkkxPxA(1.4.3)为了便于计算,在数值积分理论中已计算出各阶勒让德多项式的根及求积系数,列于下表中。高斯—勒让德求积公式的节点和系数表nkxkAnkxkA102±0.66120938650.36076157302±0.57735026921±0.23861918610.49761393463±0.77459666920.55555555567±0.94910791230.129484966200.8888888889±0.74153118560.27970539154±0.86113631160.3478548481±0.40584515140.3818300505±0.33998104360.652145154900.41795918375±0.90617984590.23692688518±0.96028985650.1012285363±0.53846931060.4786286705±0.79666647740.222381034500.5688888899±0.52553240990.31370664596±0.93246951420.1713244924±0.18343464250.3626837834对于一般的积分badxxf)(只须作代换tabbax22(1.4.4)当bxa时有11t。于是有11)22(2)(dttabbafabxdxxfba(1.4.5)这样，对上式右端的积分就可以应用高斯—勒让德求积公式了。高斯—勒让德求积公式可以推广到二维、三维的情形中，为此，把二重积分，再用一维高斯—勒让德求积公式(1.4.1)即有11111111)()(dxyxfdydxdyyxfniiiiniidyyxfAdyyxfA111111)()((ninjjijininjjijiyxAAyxfAA1111)()((于是得二维高斯—勒让德求积公式如下ninjjijiyxAAdxdyyxf111111)()((1.4.6)对于一般形式ninjjijiyxAAcdacdxdyyxf111111)(22)((1.4.7)其中22abxabxii22dcycdyij类似地还有三维高斯—勒让德求积公式如下dzdxdyzyxf111111)(ninjnkkjikjizyxAAA111)((1.4.8)这里的求积系数kjiAAA、、等仍如前所述。1.5第二类修正贝赛尔函数一、一般表达式n阶第二类修正贝赛尔函数的表达式为：10012)(1)1()2()1()1(21)(nkknnkknknkxkknxK！！！！nkxkknx2)2()]1(21)1(212[ln（1.5.1）其中，n=1,2,……,n=0时，去掉第一项有限和；是函数的对数导数，其性质如下：1、)(1)1(xxx2、101()()nrxnxxr3、1111'(1)()();((1)(1)nxxnxn　　　称为欧拉常数）。577216.0取n=0、n=1即可得)(0xK、)(1xK的表达式：20201()[(1)ln]()()22kkxxKxkk！(1.5.2)0121)2()]1()1(212[ln)1(1(1)(kkxkkxkkxxK！！(1.5.3)二、渐近表达式：当0x时：nnxnxK)2(2)1()(！（1.5.4）2ln)(0xxK(1.5.5)xxK1)(1(1.5.6)当x较大时：2211[4(21)]1(){1}2(8)kxrnkknrKxekxx！(1.5.7)2101(1)(21)1(){1}2(8)kkxrkkrKxekxx！(1.5.8)2111[4(21)]1(){1}2(8)kxrkkrKxekxx！(1.5.9)functionak0(x)doubleprecisiont,b1,y,ak,fif(x.gt.2.0)goto337t=dble(x)/3.75t=t*tb1=1.0d0+t*(3.5156229+t*(3.0899424+t*(1.2067492+t*(0.2659732*+t*(0.0360768+t*0.0045813)))))t=0.5d0*dble(x)y=t*tak=-dlog(t)*b1-0.57721566+y*(0.4227842+y*(0.23069756+y*(0.0348859*+y*(0.00262698+y*(0.0001075+y*0.0000074)))))goto338337if(x.gt.55.)thenak=0.0goto338endift=2.0d0/dble(x)f=dexp(dble(-x))/dsqrt(dble(x))ak=f*(1.25331414+t*(-0.07832358+t*(0.02189568+t*(-0.01062446*+t*(0.00587872+t*(-0.