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

【MATLAB算例】基于3节点三角形单元的矩形薄板分析 将此结构按三角形单元划分成432个三角形（X方向分成18段，Y方向分成12段），总共分成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;endXY=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);endDof=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);enddisp(前处理完成);前处理完成 单元刚度矩阵 E=2*1011;u=0.3; 平面应力问题D=E/(1-u2)*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; Bei,1 = B/2/area ; kei,1=Bei,1*D*Bei,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)+keie,1(n1,n2); end endend 单元等效节点荷载 y=(0:LY); P=(107/0.03)*y-107; 左右受变化的三角形荷载,在如图的坐标系下Re=sparse(ne,6);for i=1:ne; switch i case num2cell(1:ne/2-nx2) Pe=0 0 0 0 0 0; case num2cell(ne/2-nx2+1:ne/2) Pe=-LX*P*0,0,0,1,0,1/nx/2; case num2cell(ne/2+1:ne-nx2) Pe=0 0 0 0 0 0; otherwise Pe=-LX*P*0,0,0,1,0,1/nx/2; end Re(i,:)=Pe; end 荷载叠加 Rr=sparse(1,2*np);for i=1:ne a=Dof(i,1); b=Dof(i,2); c=Dof(i,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 Rr(DOF(n1)= Rr(DOF(n1)+Re(i,n1); endend 生成需处理的行列 cp=1:nx+1;ctype=ones(1,length(cp);ctype(nx2+1)=2;cp_all=(cp-1)*(ny+1)+1;p_stake=zeros(1,2*length(cp);for i=1:length(cp) switch ctype(i) case 2 p_stake(2*i)=2*cp_all(i); p_stake(2*i-1)=2*cp_all(i)-1; case 1 p_stake(2*i)=2*cp_all(i); p_stake(2*i-1)=; otherwise p_stake(2*i)=; p_stake(2*i-1)=2*cp_all(i)-1; endendm,j=find(p_stake=0);p_stake(:,j)=; 处理对应的行列 KK_d=KK;KK_f=KK;KK_d(p_stake,:)=;KK_d(:,p_stake)=;KK_f(:,p_stake)=;KK_f=KK_f(p_stake,:);Rr_unkown=Rr;Rr_unkown(:,p_stake)=;RR=transpose(Rr_unkown);L,U=lu(KK_d);UU=U(LRR);Rx=KK_f*UU; 数值计算部分 UU_all=UU;for i=1:length(p_stake)UU_all=UU_all(:,1:p_stake(i)-1),0,UU_all(:,p_stake(i):end);endfor i=1:np UU_info(i,:)=UU_all(2*i-1),UU_all(2*i);end 出图部分（运行后显示） 梁尺寸及荷载图 figure;set(gcf,outerposition,get(0,ScreenSize);set(gcf,name,梁的尺寸及荷载);line(-0.015,-0.01,0.06,0.06),hold online(-0.01,-0.01,0.06,0),hold online(-0.01,-0.005,0,0),hold online(-0.005,-0.015,0,0.06),hold online(0.105,0.1,0.06,0.06),hold on line(0.1,0.1,0.06,0),hold online(0.1,0.095,0,0),hold online(0.095,0.105,0,0.06),hold onrectangle(position,0,0,0.09,0.06),hold onquiver(-0.01,0.06,-0.0055,0,LineWidth,2.0,MaxHeadSize,0.8,color,k),hold onquiver(-0.01,0.05,-0.0038,0,LineWidth,2.0,MaxHeadSize,0.8,color,k),hold onquiver(-0.01,0.04,-0.002,0,LineWidth,2.0,MaxHeadSize,0.8,color,k),hold onquiver(-0.01,0.02,0.002,0,LineWidth,2.0,MaxHeadSize,0.8,color,k),hold onquiver(-0.01,0.01,0.0038,0,LineWidth,2.0,MaxHeadSize,0.8,color,k),hold onquiver(-0.01,0,0.0055,0,LineWidth,2.0,MaxHeadSize,0.8,color,k),hold onquiver(0.1,0.06,0.0055,0,LineWidth,2.0,MaxHeadSize,0.8,color,k),hold onquiver(0.1,0.05,0.0038,0,LineWidth,2.0,MaxHeadSize,0.8,color,k),hold onquiver(0.1,0.04,0.002,0,LineWidth,2.0,MaxHeadSize,0.8,color,k),hold onquiver(0.1,0.02,-0.002,0,LineWidth,2.0,MaxHeadSize,0.8,color,k),hold onquiver(0.1,0.01,-0.0038,0,LineWidth,2.0,MaxHeadSize,0.8,color,k),hold onquiver(0.1,0,-0.0055,0,LineWidth,2.0,MaxHeadSize,0.8,color,k),hold ontext(-0.0155,0.062,1000N/cm2)text(-0.0105,-0.002,1000N/cm2)text(0.1,0.062,1000N/cm2)text(0.095,-0.002,1000N/cm2)text(0,-0.002,梁宽9cm，高6cm，厚1cm。E=21011N/m2，=0.3（此图单位：m）)axis(-0.025 0.12 -0.015 0.07),hold on;网格图 figure;set(gcf,outerposition,get(0,ScreenSize);set(gcf,name,网格划分);for i=1:ne line(unit(i,1:3),unit(i,1),unit(i,4:6),unit(i,4);end axis(-0.025 0.12 -0.015 0.07),hold on;xytext=num2str(1:np);text(XY(:,1)+LX/nx/8,XY(:,2)+LY/ny/4,xytext);print(gcf,-dbitmap,model.bmp); 位移矢量图 figure;set(gcf,outerposition,get(0,ScreenSize);set(gcf,name,结点的位移矢量场);quiver(XY(:,1),XY(:,2),-UU_info(:,1),UU_info(:,2),axis(-0.025 0.12 -0.015 0.07),hold on;for i=1:ne plot(unit(i,1:3),unit(i,1),unit(i,4:6),unit(i,4),r:);endprint(gcf,-dbitmap,位移矢量场.bmp);disp(梁顶面结点);JD=13:13:247disp(梁顶面结点位移);WY=-UU_info(13);-UU_info(26);-UU_info(39);-UU_info(52);-UU_info(65);-UU_info(78);. -UU_info(91);-UU_info(104);-UU_info(117);-UU_info(130);UU_info(143);UU_info(156);. UU_info(169);UU_info(182);UU_info(195);UU_info(208);UU_info(221);UU_info(234);UU_info(247);结点位移梁顶面结点JD = 13 26 39 52 65 78 91 104 117 130 143 156 169 182 195 208 221 23