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

1、.1【MATLAB算例】基于3节点三角形单元的矩形薄板分析 将此构造按三角形单元划分成432个三角形*方向分成18段，Y方向分成12段，总共分成19*13=247个结点的有限元模型，具体步骤详细程序如下：tic;Initial_info=0.09 0.06 18 12;disp(该程序计算的是,num2str(Initial_info(3)+1),*,num2str(Initial_info(4)+1),=,. num2str(Initial_info(3)+1)*(Initial_info(4)+1),个结点的有限元模型);L*=Initial_info(1); LY=Initial_inf

2、o(2); n*=Initial_info(3);ny=Initial_info(4);ne=2*n*ny;np=(n*+1)*(ny+1);for i=1:n*+1; j=1:ny+1; Np(i,j)=j+(i-1)*(ny+1);end 生成节点编号矩阵Npfor i=1:n*+1; j=1:ny+1; *(i,j)=(i-1)*L*/n*; YY(i,j)=(j-1)*LY/ny;end*Y=reshape(*,np,1),reshape(YY,np,1);n*2=n*/2;Np1=Np(1:n*2+1,:);Np2=Np(n*2+1:end,:);for i=1:n*2*ny; if

3、 rem(i,n*2)=0 *p=n*2; yp=i/n*2; else *p=rem(i,n*2); yp=fi*(i/n*2)+1; end Dof1(i,:)=Np1(*p,yp),Np1(*p+1,yp),Np1(*p,yp+1); Dof1(i+n*2*ny,:)=Np1(*p+1,yp),Np1(*p+1,yp+1),Np1(*p,yp+1); Dof2(i,:)=Np2(*p,yp),Np2(*p+1,yp),Np2(*p+1,yp+1);Dof2(i+n*2*ny,:)=Np2(*p,yp),Np2(*p+1,yp+1),Np2(*p,yp+1);endDof=Dof1;Dof

4、2;for i=1:ne unit(i,:)=*Y(Dof(i,1),1),*Y(Dof(i,2),1),*Y(Dof(i,3),1),. *Y(Dof(i,1),2),*Y(Dof(i,2),2),*Y(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 *i=unit(i,1); yi=unit(i,4); *j=unit(i,2); yj=unit(i,5); *m=unit(i,3); ym=unit(i,6); ai=*j*ym

5、-*m*yj; aj=*m*yi-*i*ym; am=*i*yj-*j*yi; bi=yj-ym; bj=ym-yi; bm=yi-yj; ci=-(*j-*m); cj=-(*m-*i); cm=-(*i-*j); 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=

6、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-n*2) Pe=0

7、0 0 0 0 0;case num2cell(ne/2-n*2+1:ne/2) Pe=-L*P*0,0,0,1,0,1/n*/2; case num2cell(ne/2+1:ne-n*2) Pe=0 0 0 0 0 0; otherwise Pe=-L*P*0,0,0,1,0,1/n*/2; end Re(i,:)=Pe; end 荷载叠加Rr=sparse(1,2*np);for i=1:nea=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

8、; DOF(6)=2*c; for n1=1:6 Rr(DOF(n1)= Rr(DOF(n1)+Re(i,n1); endend生成需处理的行列cp=1:n*+1;ctype=ones(1,length(cp);ctype(n*2+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(

9、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);R*=KK_f

10、*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.0

11、6,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,Ma*HeadSize,0.8,colo

12、r,k),hold onquiver(-0.01,0.05,-0.0038,0,LineWidth,2.0,Ma*HeadSize,0.8,color,k),hold onquiver(-0.01,0.04,-0.002,0,LineWidth,2.0,Ma*HeadSize,0.8,color,k),hold onquiver(-0.01,0.02,0.002,0,LineWidth,2.0,Ma*HeadSize,0.8,color,k),hold onquiver(-0.01,0.01,0.0038,0,LineWidth,2.0,Ma*HeadSize,0.8,color,k),hol

13、d onquiver(-0.01,0,0.0055,0,LineWidth,2.0,Ma*HeadSize,0.8,color,k),hold onquiver(0.1,0.06,0.0055,0,LineWidth,2.0,Ma*HeadSize,0.8,color,k),hold onquiver(0.1,0.05,0.0038,0,LineWidth,2.0,Ma*HeadSize,0.8,color,k),hold onquiver(0.1,0.04,0.002,0,LineWidth,2.0,Ma*HeadSize,0.8,color,k),hold onquiver(0.1,0.0

14、2,-0.002,0,LineWidth,2.0,Ma*HeadSize,0.8,color,k),hold onquiver(0.1,0.01,-0.0038,0,LineWidth,2.0,Ma*HeadSize,0.8,color,k),hold onquiver(0.1,0,-0.0055,0,LineWidth,2.0,Ma*HeadSize,0.8,color,k),hold onte*t(-0.0155,0.062,1000N/cm2)te*t(-0.0105,-0.002,1000N/cm2)te*t(0.1,0.062,1000N/cm2)te*t(0.095,-0.002,

15、1000N/cm2)te*t(0,-0.002,梁宽9cm，高6cm，厚1cm。E=21011N/m2，=0.3此图单位：m)a*is(-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 a*is(-0.025 0.12 -0.015 0.07),hold on;*yte*t=num2str(1:np);te*t(*Y

16、(:,1)+L*/n*/8,*Y(:,2)+LY/ny/4,*yte*t);print(gcf,-dbitmap,model.bmp);位移矢量图figure;set(gcf,outerposition,get(0,ScreenSize);set(gcf,name,结点的位移矢量场);quiver(*Y(:,1),*Y(:,2),-UU_info(:,1),UU_info(:,2),a*is(-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:);endpr

17、int(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 234