【MATLAB算例】基于3節(jié)點(diǎn)三角形單元的矩形薄板分析匯總

[MATLAB算例】基于3節(jié)點(diǎn)三角形單元的矩形薄板分析將此結(jié)構(gòu)按三角形單元?jiǎng)澐殖?32個(gè)三角形(X方向分成18段，Y方向分成12段)，總共分成19X13=247個(gè)結(jié)點(diǎn)的有限元模型，具體步驟詳細(xì)程序如下：tic;Initial_info=[0.090.061812];disp(該程序計(jì)算的是',num2str(Initial_info⑶+1),'X',num2str(Initial_info(4)+1),'=',...num2str((Initial_info(3)+1)*(Initial_info(4)+1)),'個(gè)結(jié)點(diǎn)的有限元模型']);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);fori=1:nx+1;j=1:ny+1;Np(i,j)=j+(i-1)*(ny+1);end生成節(jié)點(diǎn)編號(hào)矩陣Npfori=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,:);fori=1:nx2*ny;ifrem(i,nx2)==0xp=nx2;yp=i/nx2;elsexp=rem(i,nx2);yp=fix(i/nx2)+1;endDof1(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];fori=1:neunit(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('前處理完成');前處理完成單元?jiǎng)偠染仃嘐=2*10A11;u=0.3;平面應(yīng)力問題D=E/(1-uA2)*[1u0;u10;00(1-u)/2];fori=1:nexi=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=[bi0bj0bm00ci0cj0cmcibicjbjcmbm];Be{i,1}=B/2/area;ke{i,1}=[Be{i,1}]'*D*Be{i,1}*area;end總剛度矩陣疊加KK=sparse(2*np,2*np);forie=1:nea=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;forn1=1:6forn2=1:6KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+ke{ie,1}(n1,n2);endendend單元等效節(jié)點(diǎn)荷載y=(0：LY);P=(10A7/0.03)*y-10A7;左右受變化的三角形荷載，在如圖的坐標(biāo)系下Re=sparse(ne,6);fori=1:ne;switchicasenum2cell(1:ne/2-nx2)Pe=[000000];casenum2cell(ne/2-nx2+1:ne/2)Pe=-LX*P*[0,0,0,1,0,1]/nx/2;casenum2cell(ne/2+1:ne-nx2)Pe=[000000];otherwisePe=-LX*P*[0,0,0,1,0,1]/nx/2;endRe(i,:)=Pe;end荷載疊加Rr=sparse(1,2*np);fori=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;DOF(6)=2*c;forn1=1:6Rr(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));fori=1:length(cp)switchctype(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)=口；otherwisep_stake(2*i)=[];p_stake(2*i-1)=2*cp_all(i)-1;endend[m,j]=find(p_stake==0);p_stake(:,j)=[];處理對(duì)應(yīng)的行列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\(L\RR);Rx=KK_f*UU;數(shù)值計(jì)算部分UU_all=UU';fori=1:length(p_stake)UU_all=[UU_all(:,1:p_stake(i)-1),0,UU_all(:,p_stake(i):end)];endfori=1:npUU_info(i,:)=[UU_all(2*i-1),UU_all(2*i)];end出圖部分(運(yùn)行后顯示)梁尺寸及荷載圖figure;set(gcf,'outerposition',get(0,'ScreenSize'));set(gcf,'name',梁的尺寸及荷載');line([-0.015,-0.01],[0.06,0.06]),holdonline([-0.01,-0.01],[0.06,0]),holdonline([-0.01,-0.005],[0,0]),holdonline([-0.005,-0.015],[0,0.06]),holdonline([0.105,0.1],[0.06,0.06]),holdonline([0.1,0.1],[0.06,0]),holdonline([0.1,0.095],[0,0]),holdonline([0.095,0.105],[0,0.06]),holdonrectangle('position',[0,0,0.09,0.06]),holdonquiver(-0.01,0.06,-0.0055,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),holdonquiver(-0.01,0.05,-0.0038,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),holdonquiver(-0.01,0.04,-0.002,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),holdonquiver(-0.01,0.02,0.002,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),holdonquiver(-0.01,0.01,0.0038,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),holdonquiver(-0.01,0,0.0055,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),holdonquiver(0.1,0.06,0.0055,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),holdonquiver(0.1,0.05,0.0038,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),holdonquiver(0.1,0.04,0.002,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),holdonquiver(0.1,0.02,-0.002,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),holdonquiver(0.1,0.01,-0.0038,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),holdonquiver(0.1,0,-0.0055,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),holdontext(-0.0155,0.062,'1000N/cmA2')text(-0.0105,-0.002,'1000N/cmA2')text(0.1,0.062,'1000N/cmA2')text(0.095,-0.002,'1000N/cmA2')text(0,-0.002，'梁寬9cm,高6cm,厚1cm。E=2X10A11N/mA2,仙=0.3(此圖單位:m)')axis([-0.0250.12-0.0150.07]),holdon;網(wǎng)格圖figure;set(gcf,'outerposition',get(0,'ScreenSize'));set(gcf,'name',網(wǎng)格劃分');fori=1:neline([unit(i,1:3),unit(i,1)],[unit(i,4:6),unit(i,4)]);endaxis([-0.0250.12-0.0150.07]),holdon;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',結(jié)點(diǎn)的位移矢量場');quiver(XY(:,1),XY(:,2),-UU_info(:,1),UUJnfo(:,2)),axis([-0.0250.12-0.0150.07]),holdon;fori=1:neplot([unit(i,1:3),unit(i,1)],[unit(i,4:6),unit(i,4)],'r:');endprint(gcf,'-dbitmap','位移矢量場.bmp');disp('梁頂面結(jié)點(diǎn)');JD=13:13:247disp('梁頂面結(jié)點(diǎn)位移');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);]結(jié)點(diǎn)位移梁頂面結(jié)點(diǎn)JD=1326395265789110411713014315616