MATLAB算例基于3节点三角形单元的矩形薄板分析.docx
- 文档编号:11337080
- 上传时间:2023-02-27
- 格式:DOCX
- 页数:12
- 大小:16.46KB
MATLAB算例基于3节点三角形单元的矩形薄板分析.docx
《MATLAB算例基于3节点三角形单元的矩形薄板分析.docx》由会员分享,可在线阅读,更多相关《MATLAB算例基于3节点三角形单元的矩形薄板分析.docx(12页珍藏版)》请在冰豆网上搜索。
MATLAB算例基于3节点三角形单元的矩形薄板分析
【MATLAB算例】基于3节点三角形单元的矩形薄板分析
将此结构按三角形单元划分成432个三角形(X方向分成18段,Y方向分成12段),总共分成19X13=247个结点的有限元模型,具体步骤详细程序如下:
tic;
Initial_info=[0.090.061812];
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);
fori=1:
nx+1;
j=1:
ny+1;
Np(i,j)=j+(i-1)*(ny+1);
end
生成节点编号矩阵Np
fori=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,:
);
fori=1:
nx2*ny;
ifrem(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];
fori=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)*[1u0;u10;00(1-u)/2];
fori=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=[bi0bj0bm0
0ci0cj0cm
cibicjbjcmbm];
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:
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;
forn1=1:
6
forn2=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);
fori=1:
ne;
switchi
casenum2cell(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];
otherwise
Pe=-LX*P*[0,0,0,1,0,1]/nx/2;
end
Re(i,:
)=Pe;
end
荷载叠加
Rr=sparse(1,2*np);
fori=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;
forn1=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));
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)=[];
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';
fori=1:
length(p_stake)
UU_all=[UU_all(:
1:
p_stake(i)-1),0,UU_all(:
p_stake(i):
end)];end
fori=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]),holdon
line([-0.01,-0.01],[0.06,0]),holdon
line([-0.01,-0.005],[0,0]),holdon
line([-0.005,-0.015],[0,0.06]),holdon
line([0.105,0.1],[0.06,0.06]),holdon
line([0.1,0.1],[0.06,0]),holdon
line([0.1,0.095],[0,0]),holdon
line([0.095,0.105],[0,0.06]),holdon
rectangle('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/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.0250.12-0.0150.07]),holdon;
网格图
figure;
set(gcf,'outerposition',get(0,'ScreenSize'));
set(gcf,'name','网格划分');
fori=1:
ne
line([unit(i,1:
3),unit(i,1)],[unit(i,4:
6),unit(i,4)]);
end
axis([-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','结点的位移矢量场');
quiver(XY(:
1),XY(:
2),-UU_info(:
1),UU_info(:
2)),
axis([-0.0250.12-0.0150.07]),holdon;
fori=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(234);UU_info(247);]
结点位移
梁顶面结点
JD=
13263952657891104117130143
156169182195208221234247
梁顶面结点位移
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),gridon
toc;
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- MATLAB 基于 节点 三角形 单元 矩形 薄板 分析