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

【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;

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);

P=(10^7/0.03)*y-10^7; 左右受变化的三角形荷载,在如图的坐标系下

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);

end

end

生成需处理的行列

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;

end

end

[m,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\(L\RR);

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)];

end

for 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 on

line([-0.01,-0.01],[0.06,0]),hold on

line([-0.01,-0.005],[0,0]),hold on

line([-0.005,-0.015],[0,0.06]),hold on

line([0.105,0.1],[0.06,0.06]),hold on

line([0.1,0.1],[0.06,0]),hold on

line([0.1,0.095],[0,0]),hold on

line([0.095,0.105],[0,0.06]),hold on

rectangle('position',[0,0,0.09,0.06]),hold on

quiver(-0.01,0.06,-0.0055,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(-0.01,0.05,-0.0038,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(-0.01,0.04,-0.002,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(-0.01,0.02,0.002,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(-0.01,0.01,0.0038,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(-0.01,0,0.0055,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(0.1,0.06,0.0055,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(0.1,0.05,0.0038,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(0.1,0.04,0.002,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on

quiver(0.1,0.02,-0.002,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(0.1,0.01,-0.0038,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on quiver(0.1,0,-0.0055,0,'LineWidth',2.0,'MaxHeadSize',0.8,'color','k'),hold on

text(-0.0155,0.062,'1000N/cm^2')

text(-0.0105,-0.002,'1000N/cm^2')

text(0.1,0.062,'1000N/cm^2')

text(0.095,-0.002,'1000N/cm^2')

text(0,-0.002,'梁宽9cm，高6cm，厚1cm。E=2×10^11N/m^2，μ=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:');

end

print(gcf,'-dbitmap','位移矢量场.bmp');

disp('梁顶面结点');

JD=13:13:247

disp('梁顶面结点位移');

WY=[-UU_info(13);-UU_info(26);-UU_info(39);-UU_info(52);-UU_info(65);-UU_i nfo(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(2 34);UU_info(247);]

结点位移

梁顶面结点

JD =

13 26 39 52 65 78 91 104 117 130 143

156 169 182 195 208 221 234 247 梁顶面结点位移

WY =

1.0e-006 *

0.6750

0.6000

0.5250

0.4500

0.3750

0.3000

0.2250

0.1500

0.0750

0.0000

-0.0750

-0.1500

-0.2250

-0.3000

-0.3750

-0.4500

-0.5250

-0.6000

-0.6750

figure;

set(gcf,'outerposition',get(0,'ScreenSize'));

set(gcf,'name','结点位移');

plot(JD,WY,'r:','LineWidth',2.0),grid on

toc;