弹性动力学的边界无单元法

中国科学G辑物理学力学天文学2005,35(4):435-448435弹性动力学的边界无单元法程玉民彭妙娟上海大学上海市应用数学和力学研究所,上海200072)上海大学土木工程系,上海200072)摘要讨论了 Hilbert空间上的改进的移动最小二乘法,并将弹性动力学的边界积分方程方法与改进的移动最小二乘法结合,提出了弹性动力学的边界无单元法.改进的移动最小二乘法避免了求解病态方程组,既提高了效率,又提高了精度.边界无单元法是边界积分方程的无网格方法的直接列式法,容易引入边界条件,且具有更高的精度.最后给出了弹性动力学的边界无单元法的数值实现方法和具体的算例关键词移动最小二乘法边界积分方程无网格方法弹性动力学改进的移动最小二乘法边界无单元法无网格方法近年来得到了较大发展,是目前计算力学研究的热点之无网格方法采用基于点的近似,在处理诸如大变形、裂纹扩展等问题时不需进行网格重构,具有较为明显的优越性.目前发展的无网格方法有扩散单元法( Diffuse element method,即DEM)、无单元 Galerkin法( lement-Free GalerkinMethod,即EFG)、Hp- -clouds无网格方法、有限点法( The Finite Point method,即FPM)、无网格局部 Petrov-Galerkin方法( Meshless local petrov- Galerkin Method,即MLPG)、多尺度重构核粒子方法(Mult- Scale Reproducing Kernel ParticleMethod,即MRKP)、小波粒子方法( Wavelet Particle Method,即WPM)、径向基函数法( Radial basis Functions,即RBF)、无网格有限元法( Meshless finite elementMethod,即MFEM)、移动粒子有限元法( Moving particle Finite Element Method,即 MPFEM)以及边界积分方程的无网格方法等边界积分方程的无网格方法是将边界积分方程方法与无网格方法中的移动最小二乘法结合而形成的.目前这一方面的研究主要有 Mukherjee等人提出的边界点法( Boundary Node method,即BNM)3、Atui等人提出的局部边界积分方004-08-02收稿,2005-04-22收修改稿中国煤化工*上海市重点学科建设项目资助(批准号:Y0103)*sE-mail:yecheng@sh163net:yecheng@staff.shu.edu.cnCNMHGSCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy436中国科学G辑物理学力学天文学第35卷程方法 Local Boundary Integral Equation Method,即LBI)、姚振汉等人提出的杂交边界点法( Hybrid Boundary Node Method)9、以及本文作者提出的边界无单元法( Boundary Element-Free Method,即BEFM)等与边界点法和局部边界积分方程方法相比,边界无单元法在形成的边界积分方程中采用节点变量的真实解为未知量,而边界点法和局部边界积分方程方法均采用节点变量的近似解为未知量.边界无单元法是边界积分方程无网格方法的直接解法.边界无单元法不仅直观而且具有较高的精度,还可方便地引入边界条件本文在 Hilbert空间上讨论了改进的移动最小二乘法,并将弹性动力学的边界积分方程方法与改进的移动最小二乘法结合,提出了弹性动力学的边界无单元法,最后给出了弹性动力学的边界无单元法的数值实现方法和具体的算例1移动最小二乘法在移动最小二乘法( Moving least- square approximation)中,取试函数n(x)=∑n(x)a(x)=p(x)a(x)(x∈9)(1)为函数u(x)的逼近函数.这里m是基函数的个数,P(x)是基函数,a1(x)是相应的系数.基函数的通常形式为线性基:p1=(1,x)(对一维区域)p=(,x,x2)(对二维区域二次基p2=(1,x1,x2)(对一维区域p=(1,x1,x2,x,xx2,x2)(对二维区域)对应于(1)式的整体逼近,在点x的邻域内的局部逼近定义为n'(x,)=∑(x)a1(x)=p(x)a(x),(6)其中系数a1(x)根据加权最小二乘法来确定,它使得对函数n(x)的局部逼近误差最小定义J=2w(r-x,u(r,x)-u(xn)w(x-x1)∑P(x7)a中国煤化工其中wx-x1)是具有紧支集特性的权函数,xCNMH点SCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy第4期程玉民等:弹性动力学的边界无单元法437(7)式可用矩阵形式表示为J=(Pa-u) w(r(pa-u)(8)其中(r,) p2(n)l pm(n)P=|()P2()LPm(与)(10)p,(An) P2(r,n)L Pm(n)L为了得到a(x),对J取极值,即得A(x)a(x)-B(r)u=0(12)A(x)a(r)=b(x )u(13)其中矩阵A(x)和B(x)分别为A(r)=P w(r)P(14)B(x)=Pw(x)(15)由(13)式,可得a(r)=A (xB(x)u(16)这样,逼近函数u(x)的表达式为x)=φ(x)m=∑(17)其中Φ(x)为形函数Φ(x)=闽(x,2(x),,2(x)=p(x)A(x)B(x)(18)2改进的移动最小二乘法在移动最小二乘法中,当m较大时,方程(13)有时是病态的,甚至是奇异的这样,方程(13就难以求解或获得正确的解.若选取正交函数作为基函数,则所得到的方程即不病态也不奇异,而且不需求中国煤化工方程组的解THCNMHG中国科学G辑物理学力学天文学第35卷21 Hilbert空间span(p)为了讨论正交基函数,我们先来证明pam(p)是一个 Hilbert空对vf(x),g(x)∈spmn(p),定义I=1那么(f,g)是一个内积,spmn(p)是一个内积空间.下面我们给出证明(a)由(19)式可得(f,g)=(g,f)(b)和β为常数,那么(f+βg,h)=(f,h)+β(g,h),(21)(c)对wf(x)∈spmn(p),可得1)f(x1)2≥0(22)若(f,f)=0,那么w(x-x1)=0或f(x1)=0,I=1,2,…,N.在移动最小二乘法中(x-x1)≠0、所以我们有f(x1)=0,I=1,2,…,N所以(f,g)是一个内积由于spm(p)是一个线性空间,则spmn(p)是一个内积空间对内积空间 spani(p)来说,如果它是完备的,那么它是一个 Hilbert空间.现在我们来证明spn(p)是完备的设{(x)为span(p)中的任一收敛级数,且lim fn (x)=f(r)(23)将fn(x)表示为fn(r)1,2,L其中an为常数.那么mf(x)=lm∑ankP(x)=∑(imak)(x)由于极限limf1(x)存在,则极限 lim ank存在,即有lim可得中国煤化工CNMHGSCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy第4期程玉民等:弹性动力学的边界无单元法439lim fn(xf(x)=∑akPk(x)(28)因此f(x)∈spon(p).这就证明了spn(p)是完备的,所以span(p)是一个 Hilbert空间2.2带权的正交函数族对于点集{x;}和权函数{w},若一组函数q1(x),q2(x,L,qn(x)满足如下条件0 k(k,9)=∑呢92(x)(x)(k,j=1,2,L,m)则称φ1(x),q2(x),L,φn(x)是关于点集{x;}带权{w}的正交函数族.当q(x),q2(x,L,φn(x)是多项式时,就称q1(x),q2(x),L,pn(x)是关于点集{x;}带权{w}的正交多项式2.3改进的移动最小二乘法方程(13)可写成(P1,P1)(P1,P2)L(1,Pmn)a(x)|(D1,4)(P2,P1)(P2P2)L(P2Pm)ax2(x)_(p2,a)(30)(Pm, P1)(Pm, P2)L (Pm, Pm)Lam(x)L(pm, u,)若{P(x),i=1,2,L,m,为 Hilbert空间span(p)上的关于点集{x}的带权的正交基函数族,即Pi, Pj(≠j)则方程(30)可写成0(P2,P20p2,ur(32)这样我们可以直接得到系数a(x),即.1TYH魏中国科学G辑物理学力学天文学第35卷写成矩阵形式a(x)=A(rb(x)u(34)其中(P1,P1)(P2,P2)AO(35)将(33)代入(1)式可得)=(x)=∑(x(36)其中Φ(x)为形函数Φ(x)=应1(x,2(x),,(x)=P(x)A(x)B(x)(37)这样,系数a(x)可以简单、直接地得到,不需要求矩阵A(x)的逆,避免了求解病态或奇异的方程组,既提高了效率,又提高了精度.2.4离散点上带权的正交基函数的构造利用 Schmidt正交化方法,带权的正交基函数可构造如下P11,PPki=2.3k= (pk,Pk)也可表示为P1p239)P=(r-a1)P-1-bP1其中rP;-1,P2-1)(P1-1,Pb=(41)中国煤化工对一维问题,r=x;对二维问题,r=f(YHCNMHG+x2或SCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy第4期程玉民等:弹性动力学的边界无单元法r=x+x2等;对三维问题,r=f(x,x2,x,例如可取r=√x2+x2+鸡或r=x+寺另外,也可利用 Schmidt正交化方法将类似于(2)-(5)式的基函数正交化.例如对基函数q=(q1)=(1,x,,,x2,x2…)(42)正交化得到的正交基函数为P;=q;Pk(i=1,2,3,L)若选取(38)或(39)式作为基函数,则在在试函数阶次相同的情况下,试函数中的待定常数的数目比原来要少.对线性基,原来的待定常数是3个,现在是2个;对二次基,原来的待定常数是6个,现在是3个.这样,对任一场点来说,其紧支域(影响域中所含的最小节点数就大大减少了,进而在整个求解域中所需选取的节点数目也可以大大减少.所以在精度相同的情况下,改进的移动最小二乘法形成的无网格方法比移动最小二乘法形成的无网格方法的节点数目要少得多3弹性动力学的边界无单元法对弹性动力学的基本方程进行 Fourier变换,得到 Fourier变换域中的基本方程.在 Fourier变换域中应用边界无单元法进行求解,然后应用数值 Fourier本征反变换即得时间域中的解弹性动力学的基本方程:运动微分方程(C2-C2)1(x,D)+Clu1(x,1)+f(x,)=(x,1)(x∈92),(44)其中+2pC(45)分别是弹性体内膨胀波和畸变波的传播速度,λ和μ是Lame常数,ρ是体密度物理方程o(x,)=p(C2-2C2)xmm(x,16+C2(n1(x,1)+m1(x,1)(x∈9),(46)其中δ为 Kronecker应力和位移满足如下的边界条件t, (r, t)=o(x,)中国煤化工(47)l(x,1)=q1(x,1)CNMHGhina co442中国科学G辑物理学力学天文学第35卷和初始条件(x∈),其中r和r分别表示已知面力和位移的边界,r∪I"=,⌒r"=Φ;n(x)为r上点x处的外法线的方向余弦;和分别为已知的位移和面力分量;u0和io分别为弹性体的初始位移和初始速度场将 Fourier变换f(x,o)=FLf(x,D1=6-f(r, t) - dr作用于时间域中的弹性动力学基本方程,并令l(x,O)=F[L1(x,1)](r, o)=Flo,(r, t)t;(x,0)=F[t1(x,t)]f (x, o)=FIf, (r, t)],Pi(FLP; (x, t)l,即得 Fourier变换域中的弹性动力学基本方程:运动微分方程(C2-C2x0(x,0)+C1(x,0)+f(x,0)=02(x,0)(x∈9),(51)边界条件1(x,0)=G(x,0)m(x)=(x,0)(x∈r(52)u(x,0)=q;(x,O)物理方程0(x,0)=p[(C1-2C2)mm(x,O)6+C2(1(x,)+:(x,0)(x∈9)(53)为方便起见,我们假设初始条件和体力均为零.这样,由加权残数法可得弹性动力学问题在 Fourier变换域中的边界积分方程lCn(q)x(q)=U(p,9)(p)dI-」T(P,q)(p)dr其中q为边界点源点,P为场点,Cn(q)为自由项,U和T为 Fourier变换域中弹性动力学方程的基本解.将边界r离散为N个边界子域,这些边于,亚数无关,仅供积分时使用,则边界积分方程(54)变为THE中国煤化工CNMHGSCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy第4期程玉民等:弹性动力学的边界无单元法C(o(=∑J1(p(p-∑∫Tk(p, q)u, (p)dr在各边界子域中配置一定数量的点,并分别选取相应的紧支域,这些紧支域可以相交,但其并集需覆盖整个边界.由改进的移动最小二乘法得插值公式∑Φ(x∑(x)其中u;l=u(r)L r=t;(x将插值公式(56和(57代入(55)式,则有C()(q)=∑∫(,9∑画(p)qmR(P,41)∑(P)(q(60)其中q为节点,n1个q的紧支域覆盖p,所选用的权函数是Gaus函数(d≤1)0其中c是常数,d是紧支域大小的度量.通常对二维问题的边界区域,紧支域为一维区域;对三维问题的边界区域,紧支域可选规则的圆域或矩形域对离散化的边界积分方程(60)进行数值积分后,写成矩阵形式即为([C]+[H])U]=[G][T其中U]=[1,12,l13,21,21T中国煤化x(660将边界条件代入后,对线性代数方程组(HYH中边界节CNMHG44中国科学G辑物理学力学天文学第35卷点的位移和面力变换域中内点x的位移可由内点的离散的边界积分方程∑∫rU(pq∑画(p)∑nr(p9,(p()d(67)得到,然后由几何方程和物理方程即可得到变换域中内点的应力.在求得变换域中的解后,利用数值 Fourier本征反变换即可得到时间域中的解.这就是弹性动力学的边界无单元法4数值 Fourier本征反变换ll与 Fourier变换(49)式对应的反变换为f(x,)=FI(,0)=6_f(x,o)e'o'do(68)由 Fourier本征变换2可知,(68)式可以表示为f(x,1)=∑7a1甲(),其中(7(x,0)9n「f(x.)n(o)doO)qn()为本征基函数9()=/~1(71)由(68)式可得数值 Fourier本征反变换的表达式为(2r2f(x,1)=∑"anh1(e2=∑anhn()e2(72)其中dh, (t)∑4(f(x,0)中国煤化工CNMHGSCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy第4期程玉民等:弹性动力学的边界无单元法5弹性动力学边界无单元法的数值实现对上述导出的弹性动力学边界无单元法,其数值实现过程如下(i)在问题所在域的边界配置n个节点(ⅱi)确定各节点的紧支域rn的大小d,使得UIn=In(i)选取基函数和权函数;(iv)计算形函数(ⅴ)将边界离散为N个积分子域In,得到离散的边界积分方程(ⅵi)通过数值积分,并代入边界条件,得到线性代数方程组;(ⅶi)求解此线性代数方程组,即得各边界节点的位移和面力(ⅷⅲ)由改进的移动最小二乘法的插值公式得到其他边界点的位移和面力(ⅸx)由内点的离散的边界积分方程计算域内任意一点的位移(x)由几何方程和物理方程计算域内任意一点的应力(xi)由数值 Fourier本征反变换求时间域中的解6数值算例以下采用本文提出的弹性动力学的边界无单元法对受突加荷载作用的中心圆孔板和中心裂纹板进行计算61中心圆孔板长为36cm、宽为20cm的矩形板,中央有一直径是10cm的圆形孔洞,在平面应力情况下两侧当t=0时作用突加荷载p=75MNm2,如图1所示.材料的弹性模量E=2.1×105MNm2, Poisson比=0.3,质量密度p=0785kNm2x2中国煤化工图1受突加荷载作用哨YHCNMHGhina. com446中国科学G辑物理学力学天文学第35卷我们用边界无单元法对本例进行了计算.由于对称性,我们仅考虑此矩形裂纹板的四分之一情形.边界无单元法的边界节点设置如图2所示计算所得的点A的位移值与有限元法计算结果的比较如图3所示.可以看出,本文讨论的边界无单元法的计算结果和图2中心圆孔板的边界节点配置有限元法的计算结果吻合较好.0.00400020.008边界无单元法0.012有限元法0.0010.0020.0030.004图3点A处的位移62中心裂纹板受突加荷载作用的矩形裂纹板如图4所示.裂纹长度为2a,a=0.24cm.材料的弹性模量E=2.0×103MNm2, Poisson比v=0.3,质量密度p=0.005MN.s2/n我们用边界无单元法对本例进行了计算.由于对称性,我们仅考虑此矩形裂纹板的四分之一情形.边界无单元法的边界节点设置如图5所示计算所得的正则动态应力强度因子(K()/K1,K1为静态应力强度因子)与边界元法的计算结果的比较如图6所示由图6可以看出,边界无单元法的计算结果与边界元法的计算结果吻合得很好中国煤化工CNMHGSCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy第4期程玉民等:弹性动力学的边界无单元法图4中心裂纹板图5中心裂纹板的边界节点配置3.02.50.50.0边界元单元法图6正则应力强度因子7结论和讨论本文讨论了 Hilbert空间上的改进的移动最小二乘法.改进的移动最小二乘法不会形成奇异或病态的方程组,可提高计算精度和效率改进的移动最小二乘法中所含的待定常数较少,这样只需较少的节点的变量的值就可通过逼近函数得到任意场点的变量的值,那么由改进的移动最小二乘法形成的无网格方法求解问题时可大大减少中国煤化工本文建立了求解弹性动力学问题的边界eHCNMHG公式,提hina co中国科学G辑物理学力学天文学第35卷出了具体的数值实现方法.数值算例表明,本文方法是有效的本文方法可推广到其他边界元法可以解决的问题.参考文献I Belytschko T, Krongauz Y, Organ D, et al. Meshless methods: An overview and recent developments.Computer Methods in Applied Mechanics and Engineering, 1996. 139: 3-472 Li S F, Liu W K. Meshfree and particle methods and their applications. Applied Mechanics Review2002,55:1-343 Mukherjee Y X, Mukherjee S The boundary node method for potential problems. International Journal foNumerical Methods in Engineering, 1997, 40: 797-8154 Kothnur V S, Mukherjee S, Mukherjee Y X. Two dimensional linear elasticity by the boundary nodemethod, International Journal of Solids and Structures, 1999: 36: 1129~1145 Chati M K, Mukherjee S, Mukherjee Y x. The boundary node method for three-dimensional linearelasticity, International Journal for Numerical Methods in Engineering, 1999. 46: 1163-11846 Chati M K, Mukherjee S. The boundary node method for three-dilal problems in potential theorInternational Journal for Numerical Methods in Engineering, 2000, 47: 1523-1547 Zhu T L, Zhang J D, Atluri S N, A meshless numerical method based on the local boundary integralequation (LBIE)to solve linear and non-linear boundary value problems. Engineering Analysis withBoundary Elements, 1999, 23: 375-3898 Atluri S N, Sladek J, Sladek V, et al. The local boundary integral equation(LBIE)and it's meshlessiplementation for linear elasticity. Computational Mechanics, 2000, 25: 180-199 Zhang J M, Yao Z H, Li H. A hybrid boundary node method. International Journal for NumericalMethods in Engineering, 2002, 53: 751-76310程玉民,陈美娟.弹性力学的一种边界无单元法.力学学报,2003,35(2):181-1861嵇醒,臧跃龙,程玉民.边界元法进展及通用程序.上海:同济大学出版社,199712陈绍汀,韩海潮.傳立叶本征变换FI.应用力学学报,1987,4(1):33~38中国煤化工CNMHGSCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy