2Dmtlab.rar_2Dmtlab_meshless_悬臂梁_悬臂梁无网格_梁

clear; Lb=48; D=12; young=30e6;nu=0.3; P=1000; Dmat=(young/(1-nu^2))*[1 nu 0;nu 1 0;0 0 (1-nu)/2]; ndivl=10;ndivw=4; [x,conn1,conn2,numnod]=mesh2(Lb,D,ndivl,ndivw); dmax=3.5; xspac=Lb/ndivl; yspac=D/ndivw; dm(1,1:numnod)=dmax*xspac*ones(1,numnod); dm(2,1:numnod)=dmax*yspac*ones(1,numnod); ndivlq=10;ndivwq=4; [xc,conn,numcell,numq]=mesh2(Lb,D,ndivlq,ndivwq); quado=4; [gauss]=gauss2(quado); numq2=numcell*quado^2; gs=zeros(4,numq2); [gs]=egauss(xc,conn,gauss,numcell); k=sparse(numnod*2,numnod*2); for gg=gs gpos=gg(1:2); weight=gg(3); jac=gg(4); 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); for j=1:L Bmat(1:3,(2*j-1):2*j)=[dphix(j) 0;0 dphiy(j);dphiy(j) dphix(j)]; end for i=1:L en(2*i-1)=2*v(i)-1; en(2*i)=2*v(i); end k(en,en)=k(en,en)+sparse((weight*jac)*Bmat'*Dmat*Bmat); end ind1=0;ind2=0; for j=1:numnod if(x(1,j)==0.0) ind1=ind1+1; nnu(1,ind1)=x(1,j); nnu(2,ind1)=x(2,j); end if(x(1,j)==Lb) ind2=ind2+1; nt(1,ind2)=x(1,j); nt(2,ind2)=x(2,j); end end lthu=length(nnu); ltht=length(nt); ubar=zeros(lthu*2,1); f=zeros(numnod*2,1); ind=0; gauss=gauss2(quado); for i=1:(ltht-1) ycen=(nt(2,i+1)+nt(2,i))/2; jcob=abs((nt(2,i+1)-nt(2,i))/2); for j=1:quado mark(j)=ycen-gauss(1,j)*jcob; ind=ind+1; gst(1,ind)=nt(1,i); gst(2,ind)=mark(j); gst(3,ind)=gauss(2,j); gst(4,ind)=jcob; end end gsu=gst; gsu(1,1:ind)=zeros(1,ind); qk=zeros(1,2*lthu); Imo=(1/12)*D^3; for gt=gst gpos=gt(1:2); weight=gt(3); jac=gt(4); v=domain(gpos,x,dm,numnod); L=length(v); en=zeros(1,2*L); force=zeros(1,2*L); [phi,dphix,dphiy]=shape(gpos,dmax,x,v,dm); tx=0; ty=-(P/(2*Imo))*((D^2)/4-gpos(2,1)^2); for i=1:L en(2*i-1)=2*v(i)-1; en(2*i)=2*v(i); force(2*i-1)=tx*phi(i); force(2*i)=ty*phi(i); end f(en)=f(en)+jac*weight*force'; end GG=zeros(numnod*2,lthu*2); ind1=0;ind2=0; for i=1:(lthu-1) ind1=ind1+1; m1=ind1;m2=m1+1; y1=nnu(2,m1);y2=nnu(2,m2); len=y1-y2; for j=1:quado ind2=ind2+1; gpos=gsu(1:2,ind2); weight=gsu(3,j); jac=gsu(4,j); fac11=(-P*nnu(2,m1))/(6*young*Imo); fac2=P/(6*young*Imo); xp1=gpos(1,1); yp1=gpos(2,1); uxex1=(6*Lb-3*xp1)*xp1+(2+nu)*(yp1^2-(D/2)^2); uxex1=uxex1*fac11; uyex1=3*nu*yp1^2*(Lb-xp1)+0.25*(4+5*nu)*xp1*D^2+(3*Lb-xp1)*xp1^2; uyex1=uyex1*fac2; N1=(gpos(2,1)-y2)/len;N2=1-N1; qk(2*m1-1)=qk(2*m1-1)-weight*jac*N1*uxex1; qk(2*m1)=qk(2*m1)-weight*jac*N1*uyex1; qk(2*m2-1)=qk(2*m2-1)-weight*jac*N2*uxex1; qk(2*m2)=qk(2*m2)-weight*jac*N2*uyex1; v=domain(gpos,x,dm,numnod); [phi,dphix,dphiy]=shape(gpos,dmax,x,v,dm); L=length(v); for n=1:L G1=-weight*jac*phi(n)*[N1 0;0 N1]; G2=-weight*jac*phi(n)*[N2 0;0 N2]; c1=2*v(n)-1;c2=2*v(n);c3=2*m1-1;c4=2*m1; c5=2*m2-1;c6=2*m2; GG(c1:c2,c3:c4)= GG(c1:c2,c3:c4)+G1; GG(c1:c2,c5:c6)= GG(c1:c2,c5:c6)+G2; end end end f=[f;zeros(lthu*2,1)]; f(numnod*2+1:numnod*2+2*lthu,1)=-qk'; m=sparse([k GG;GG' zeros(lthu*2)]); d=m\f; u=d(1:2*numnod); for i=1:numnod u2(1,i)=u(2*i-1); u2(2,i)=u(2*i); end displ=zeros(1,2*numnod); ind=0; for gg=x ind=ind+1; gpos=gg(1:2); v=domain(gpos,x,dm,numnod); [phi,dphix,dphiy]=shape(gpos,dmax,x,v,dm); displ(2*ind-1)=phi*u2(1,v)'; displ(2*ind)=phi*u2(2,v)'; end ind=0; enorm=0; for gg=gs ind=ind+1; gpos=gg(1:2); weight=gg(3); jac=gg(4); 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); for j=1:L Bmat(1:3,(2*j-1):2*j)=[dphix(j) 0;0 dphiy(j);dphiy(j) dphix(j)]; end for i=1:L en(2*i-1)=2*v(i)-1; en(2*i)=2*v(i); end stress(1:3,ind)=Dmat*Bmat*u(en); stressex(1,ind)=(1/Imo)*P*(Lb-gpos(1,1))*gpos(2,1); stressex(2,ind)=0; stressex(3,ind)=-0.5*(P/Imo)*(0.25*D^2-gpos(2,1)^2); err=stress(1:3,ind)-stressex(1:3,ind); err2=weight*jac*(0.5*(inv(Dmat)*err)'*(err)); enorm=enorm+err2; end enorm=sqrt(enorm); displx=displ(1:2:109); disply=displ(2:2:110); displtot=[displx;disply]; function[x,conn,numcell,numq]=mesh2(length,height,ndivl,ndivw) numcell=ndivw*ndivl; numq=(ndivl+1)*(ndivw+1); for i=1:(ndivl+1) for j=1:(ndivw+1) x(1,((ndivw+1)*(i-1)+j))=(length/ndivl)*(i-1); x(2,((ndivw+1)*(i-1)+j))=-(height/ndivw)*(j-1)+height/2; end end for j=1:ndivl for i=1:ndivw elemn=(j-1)*ndivw+i; nodet(elemn,1)=elemn+(j-1); nodet(elemn,2)=nodet(elemn,1)+1; nodet(elemn,3)=nodet(elemn,2)+ndivw+1; nodet(elemn,4)=nodet(elemn,3)-1; end end conn=nodet'; function [gs]=egauss(xc,conn,gauss,numcell) index=0; one=ones(1,4); psiJ=[-1,+1,+1,-1];etaJ=[-1,-1,+1,+1]; l=size(gauss); l=l(2); for e=1:numcell for j=1:4 je=conn(j,e);xe(j)=xc(1,je);ye(j)=xc(2,je); end for i=1:l for j=1:l index=index+1; eta=gauss(1,i);psi=gauss(1,j); N=.25*(one+psi*psiJ).*(one+eta*etaJ); NJpsi=.25*psiJ.*(one+eta*etaJ); NJeta=.25*etaJ.*(one+psi*psiJ); xpsi=NJpsi*xe';ypsi=NJpsi*ye';xeta=NJeta*xe';yeta=NJeta*ye'; jcob=xpsi*yeta-xeta*ypsi; xq=N*xe';yq=N*ye'; gs(1,index)=xq; gs(2,index)=yq; gs(3,index)=gauss(2,i)*gauss(2,j); gs(4,index)=jcob; end end end function v=gauss2(k) v(1,1)=-.861136311594052575224; v(1,2)=-.339981043584856264803; v(1,3)=-v(1,2); v(1,4)=-v(1,1); v(2,1)=.347854845137453857373; v(2,2)=.652145154862546142627; v(2,3)=v(2,2); v(2,4)=v(2,1); function v=domain(gpos,x,dm,nnodes) dif=gpos*ones(1,nnodes)-x; a=dm-abs(dif); b=[a;zeros(1,nnodes)]; c=all(b>=-100*eps); v=find(c); function [phi,dphix,dphiy]=shape(gpos,dmax,x,v,dm) L=length(v); won=ones(1,L); nv=x(1:2,v); p=[won;nv]; dif=gpos*won-nv; t=dm(1:2,v)/dmax; [w,dwdx,dwdy]=cubwgt(dif,t,v,dmax,dm); B=p.*[w;w;w]; pp=zeros(3); aa=zeros(3); daax=zeros(3); daay=zeros(3); for i=1:L pp=p(1:3,i)*p(1:3,i)'; aa=aa+w(1,i)*pp; daax=daax+dwdx(1,i)*pp; daay=daay+dwdy(1,i)*pp; end pg=[1 gpos']; [L,U,PERM]=lu(aa); for i=1:3 if i==1 C=PERM*pg'; elseif i==2 C=PERM*([0 1 0]'-daax*gam(1:3,1)); elseif i==3 C=PERM*([0 0 1]'-daay*gam(1:3,1)); end D1=C(1); D2=(C(2)-L(2,1)*D1); D3=(C(3)-L(3,1)*D1-L(3,2)*D2); gam(3,i)=D3/U(3,3); gam(2,i)=(D2-U(2,3)*gam(3,i))/(U(2,2)); gam(1,i)=(D1-U(1,2)*gam(2,i)-U(1,3)*gam(3,i))/U(1,1); end phi=gam(1:3,1)'*B; dbx=p.*[dwdx;dwdx;dwdx]; dby=p.*[dwdy;dwdy;dwdy]; dphix=gam(1:3,2)'*B+gam(1:3,1)'*dbx; dphiy=gam(1:3,3)'*B+gam(1:3,1)'*dby; function [w,dwdx,dwdy]=cubwgt(dif,t,v,dmax,dm) l=length(v); for i=1:l drdx=sign(dif(1,i))/dm(1,v(i)); drdy=sign(dif(2,i))/dm(2,v(i)); rx=abs(dif(1,i))/dm(1,v(i)); ry=abs(dif(2,i))/dm(2,v(i)); if rx>0.5 wx=(4/3)-4*rx+4*rx*rx-(4/3)*rx^3; dwx=(-4+8*rx-4*rx^2)*drdx; elseif rx<=0.5 wx=(2/3)-4*rx*rx+4*rx^3; dwx=(-8*rx+12*rx^2)*drdx; end if ry>0.5 wy=(4/3)-4*ry+4*ry*ry-(4/3)*ry^3; dwy=(-4+8*ry-4*ry^2)*drdy; elseif ry<=0.5 wy=(2/3)-4*ry*ry+4*ry^3; dwy=(-8*ry+12*ry^2)*drdy; end w(i)=wx*wy; dwdx(i)=wy*dwx; dwdy(i)=wx*dwy end