MATLAB算例基于3节点三角形单元的矩形薄板分析.doc

【MATLAB算例】基于3节点三角形单元的矩形薄19X13=247个结点的有限元模型 tic; Initial_info=[0.09 0.06 18 12]; disp(该程序计算的是,num2str(Initial_info(3)+1),X,num2str(Initial_info(4)+1),=,... num2str((Initial_info(3)+1)*(Initial_info(4)+1)),个结点的有限元模型]); LX=Initial_info(1); LY=Initial_info(2); nx=Initial_info(3); ny=Initial_info(4); ne=2*nx*ny; np=(nx+1)*(ny+1); for i=1:nx+1; j=1:ny+1; Np(i,j)=j+(i-1)*(ny+1); end 生成节点编号矩阵Np for i=1:nx+1; j=1:ny+1; XX(i,j)=(i-1)*LX/nx; YY(i,j)=(j-1)*LY/ny; end XY=[reshape(XX,np,1),reshape(YY,np,1)]; nx2=nx/2; Np1=Np(1:nx2+1,:); Np2=Np(nx2+1:end,:); for i=1:nx2*ny; if rem(i,nx2)==0 xp=nx2; yp=i/nx2; else xp=rem(i,nx2); yp=fix(i/nx2)+1; end Dof1(i,:)=[Np1(xp,yp),Np1(xp+1,yp),Np1(xp,yp+1)]; Dof1(i+nx2*ny,:)=[Np1(xp+1,yp),Np1(xp+1,yp+1),Np1(xp,yp+1)]; Dof2(i,:)=[Np2(xp,yp),Np2(xp+1,yp),Np2(xp+1,yp+1)]; Dof2(i+nx2*ny,:)=[Np2(xp,yp),Np2(xp+1,yp+1),Np2(xp,yp+1)]; end Dof=[Dof1;Dof2]; for i=1:ne unit(i,:)=[XY(Dof(i,1),1),XY(Dof(i,2),1),XY(Dof(i,3),1),... XY(Dof(i,1),2),XY(Dof(i,2),2),XY(Dof(i,3),2)]; end disp(前处理完成); 前处理完成 单元刚度矩阵 E=2*10^11;u=0.3; 平面应力问题 D=E/(1-u^2)*[1 u 0;u 1 0;0 0 (1-u)/2]; for i=1:ne xi=unit(i,1); yi=unit(i,4); xj=unit(i,2); yj=unit(i,5); xm=unit(i,3); ym=unit(i,6); ai=xj*ym-xm*yj; aj=xm*yi-xi*ym; am=xi*yj-xj*yi; bi=yj-ym; bj=ym-yi; bm=yi-yj; ci=-(xj-xm); cj=-(xm-xi); cm=-(xi-xj); area=abs((ai+aj+am)/2); B = [bi 0 bj 0 bm 0 0 ci 0 cj 0 cm ci bi cj bj cm bm]; Be{i,1} = B/2/area ; ke{i,1}=[Be{i,1}]*D*Be{i,1}*area; end 总刚度矩阵叠加 KK=sparse(2*np,2*np); for ie=1:ne a=Dof(ie,1); b=Dof(ie,2); c=Dof(ie,3); DOF(1)=2*a-1; DOF(2)=2*a; DOF(3)=2*b-1; DOF(4)=2*b; DOF(5)=2*c-1; DOF(6)=2*c; for n1=1:6 for n2=1:6 KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+ke{ie,1}(n1,n2); end end end 单元等效节点荷载 y=(0:LY);