无网格数值求解方法学习小结

无网格数值求解方法——学习小结一、无网格法的介绍有限元法存在的那些问题都来源于网格，在用有限元方法处理诸如金属冲压成型、高速冲击、动态裂纹扩展、流固耦合等涉及大变形和移动边界的问题时，由于网格可能发生严重扭曲，往往需要网格重构，不但精度受到了严重影响，计算也大幅度提高，因此有限元方法在这些领域的应用遇到了困难。直接在有限元基础上对其进行改进，效果自然不会达到最好，于是研究者把革命的对象锁定在了网格上。几经尝试以后，一种基于点集的插值方法被研究者广泛采用，现今的无网格方法，一般就指的是这一类基于点集的数值方法。无网格方法的位移函数是在点的领域内构造的，并且这些区域是可以重叠的，因此在处理大变形和移动边界等问题时，没有网格的初始划分和重构问题，这不仅有利于这类问题计算精度的提高，还可以减少数值计算难度。目前已存在十余种无网格方法，它们之间的区别主要在于试函数的选择和微分方程的等效形式。虽然无网格方法对于大变形和移动边界问题具有优势，但其存在收敛性、数值稳定性和效率等问题，因此无网格方法还只能作为有限元方法的补充。无网格方法基本思想是将有限元法中的网格结构去除，完全代之以一系列的结点排列。二、求解方法方法基于位移最小二乘(MLS)近似方法—EFG(Element-freeGalerkinMethod,Belytschko,1994)。EFG方法计算稳定，精度较高，是无网格方法中较为成熟的一种方法。无网格法就目前来说，仍没有有限元法发展得那么快。而且，大规模地使用无网格法将大大增加计算时间。因此通常只需要在那些不连续、大变形或应力集中区域使用无网格法进行离散，如冲击区域、裂纹扩展区域、大变形区域等，其余区域仍然可采用其他数值方法。微分方程组0(**)TTCUDudEVFud2.1加权余量法求解域Ω中，若场函数是精确解，则在域Ω中任一点都满足微分方程，同时在边界上任一点都满足边界条件式，此时等效积分形式或等效积分弱形式必()()()0AuAuAu12=在内...12()()()0...BuBuBu在上边界条件等效积分形式0(*)TTUAudVBud等效积分弱形式然严格地得到满足。但是对于复杂的实际问题，这样的精确解往往是很难找到的，因此,人们需要设法找到具有一定精度的近似解。设u是一个近似解，即为试函数，它可以表示成为一组已知函数或Ritz基函数i的线性组合，即1nTTiiiuxaa式中ia为待定系数或Ritz基坐标。将权函数代入加权余量积分式，由于系数jb的任意性，有0,1,2,,TTTTjAjBRadRadjm上式给出了m个方程。用于求解n个待定系数ia。如果mn，则上式是超定的，需要借助于最小二乘法解。对上式进行分部积分得到等效积分弱形式的近似形式0TTTTjjCDadEFad2.2伽辽金法按照对权函数的不同选择就得到不同的加权余量的计算方法并赋以不同的名称。如果取权函数与试函数相同，则称为Galerkin方法。110nnTTTTjAiijBiiiiRadRad我们将会看到，在很多情况下，采用伽辽金法得到的求解方程的系数矩阵是对称的，这是在用加权余量法建立有限元格式时几乎毫不例外地采用伽辽金法的主要原因，而且当存在相应的泛函数时，伽辽金法与变分法往往导致同样的结果。2.3移动最小二乘近似构造方法：考虑求解域，其中共有N个结点ix(1,2,,)iN，在各个结点处有0()iiuxu，但()iiuxu。考虑计算点x（对于无网格配点法为结点；对于伽辽金无网格方法为高斯积分点），其邻域x内的近似函数可以写为1(,)()()()()mTiiiuxxpxaxpxax式中：()ipx为Rits基函数，()iax为Rits基坐标或待求系数，[]xxyz是计算点x邻域x内任意点的坐标，它包括x，m是基函数的个数。而12()[()()()]Tmpxpxpxpx，12()[()()()]Tmaxaxaxax值得注意的是，在经典Ritz方法中，Ritz基坐标是常数，并且基函数要满足位移边界条件。在式(1)中，基函数要满足如下条件：1()1()()nipxpxC式中：1,2,,im，()nC表示在域内具有直到n阶连续导数的函数空间。2.4边界条件无网格方法的结点形函数多数都不满足关系jijNix，因此位移边界条件的处理是比较困难的。若采用紧支径向基函数来构造形函数，则可以像一般有限元方法那样来处理位移边界条件。在MLS近似中，若选奇异函数为权函数，则近似函数具有插值特性即jijNix，因此可以直接施加本质边界条件。对与其他情况，可以借助拉格朗日乘子方法来处理边界条件。拉格朗日乘子法包括两种，一种是利用边界积分中直接引入边界条件，即0udddTTTTTεσufupuλ+λu-u三、具体算例左端固定的悬臂梁，右端面受抛物线剪切载荷作用主程序：ticclear;Lx=20;Ly=10;young=210;nu=0.3;q=-1;a=0;nx=30;ny=20;ndivl=10;ndivw=6;dmax=2.89;Dmat=(young/(1-nu^2))*[1nu0;nu10;00(1-nu)/2];[x,numnod,dm]=mesh1(Lx,Ly,nx,ny,dmax);figureholdonplot(x(1,1:(ny+1)),x(2,1:(ny+1)),'k-','linewidth',3);axisequal;plot(x(1,(ny+1):(ny+1):numnod),x(2,(ny+1):(ny+1):numnod),'k-','linewidth',2);plot(x(1,numnod:-1:(numnod-ny)),x(2,numnod:-1:(numnod-ny)),'k-','linewidth',2);plot(x(1,1:(ny+1):(numnod-ny)),x(2,1:(ny+1):(numnod-ny)),'k-','linewidth',2);%plot(x(1,:),x(2,:),'k.');___axisoff;plot(x(1,:),x(2,:),'k.');axisequal;axisoff;holdoff[xc,conn,numcell,numq]=mesh2(Lx,Ly,ndivl,ndivw);[nnu,nnt,numT1,numT2]=mesh3(numq,xc,Lx,Ly,a);%nnu---%nnt---%numT1--%numT2--%numq-----quado=4;[gauss]=gauss2(quado);numq2=numcell*quado^2;gs=zeros(4,numq2);[gs]=egauss(xc,conn,gauss,numcell);[k]=kjuzhen(numnod,gs,x,dm,dmax,Dmat);rfa=400e12;[ka]=kajuzhen(numnod,nnu,numT1,xc,gauss,x,dm,dmax,rfa);K=k+ka;[f]=fjuzhen(numnod,nnt,numT2,xc,gauss,x,dm,dmax,q,Ly);%fa=zeros(2*numnod,1);%[fa]=fajuzhen(nu,young,q,numnod,nnu,numT1,xc,gauss,x,dm,dmax,rfa,Ly);[fa]=fajuzhen(nu,young,q,numnod,nnu,numT1,xc,gauss,x,dm,dmax,rfa,Lx,Ly)F=f+fa;u=zeros(2*numnod,1);fori=1:numnodu2(1,i)=u(2*i-1);u2(2,i)=u(2*i);endnx1=2;ny1=10;I=Ly^3/12;fori=1:(ny1+1)xjm(1,i)=Lx/2;xjm(2,i)=-(Ly)/ny1*(i-1)+Ly;yjm(i)=-(Ly/ny1)*(i-1)+Ly/2;stress11ex(i)=-q*(Lx-xjm(1,i))*yjm(i)/I;stress12ex(i)=q/(2*I)*(Ly^2/4-yjm(i)^2);endind=0;enorm=0;forgg=xjmind=ind+1;gpos=gg(1:2);v=domain(gpos,x,dm,numnod);L=length(v);en=zeros(1,2*L);[phi,dphix,dphiy]=shape(gpos,dmax,x,v,dm);Bmat=zeros(3,2*L);forj=1:LBmat(1:3,(2*j-1):2*j)=[dphix(j)0;0dphiy(j);dphiy(j)dphix(j)];endfori=1:Len(2*i-1)=2*v(i)-1;en(2*i)=2*v(i);endstress(1:3,ind)=Dmat*Bmat*u(en);%stressex(1,ind)=;%stressex(2,ind)=0;%stressex(3,ind)=0;%err=stress(1:3,ind)-stressex(1:3,ind);%err2=weight*jac*(0.5*(inv(Dmat)*err)'*(err));%enorm=enorm+err2;end%uex=zeros(2,numnod);I=Ly^3/12;ind4=0;fori=1:numnodif(x(2,i)==Ly/2)ind4=ind4+1;uex2(ind4)=q/(6*young*I)*(3*nu*(x(2,i)-Ly/2)^2*(Lx-x(1,i))+(4+5*nu)*(Ly/2)^2*x(1,i)+(3*Lx-x(1,i))*x(1,i)^2);endfigureholdonplot(x(1,(ny+1)/2:(ny+1):numnod),u2(2,(ny+1)/2:(ny+1):numnod),'r.');plot(x(1,(ny+1)/2:(ny+1):numnod),uex2,'-');%plot(xz,u2jy,'o');xlabel('x/m','fontweight','bold');ylabel('ux/m','fontweight','bold');legend('Uynode','ExactSolution');holdoff%figure%holdon%plot(xjm(2,1:(ny1+1)),stress(1,1:ind),'r*');%plot(xjm(2,1:(ny1+1)),stress11ex(1,1:ind),'.-');%legend('EFGSolution','exactsolution');%%xlabel('y/m','fontweight','bold');%ylabel('Stress','fontweight','bold');%%holdoff%%figure%holdon%plot(xjm(2,1:(ny1+1)),stress(3,1:ind),'r*');%plot(xjm(2,1:(ny1+1)),stress12ex(1,1:ind),'.-');%legend('EFGSolution','exactsolution');%xlabel('y/m','fontweight','bold');%ylabel('Stress','fontweight','bold');holdoffToc矩形区域内均匀节点布置：解析解与无网格近似的比较：