平面三角形单元Matlab程序.doc
- 文档编号:2528486
- 上传时间:2022-10-31
- 格式:DOC
- 页数:4
- 大小:24KB
平面三角形单元Matlab程序.doc
《平面三角形单元Matlab程序.doc》由会员分享,可在线阅读,更多相关《平面三角形单元Matlab程序.doc(4页珍藏版)》请在冰豆网上搜索。
%变量说明
%NPOINNELEMNVFIXNFORCENNODE
%总结点数,单元数,受约束边界点数,结点力数,单元结点数
%COORDLNODSYOUNGPOISSTHICK
%结构结点坐标数组,单元定义数组,弹性模量,泊松比,厚度
%FORCEFIXEDBMATXDMATXSMATX
%节点力数组,约束信息数组,单元应变矩阵,单元弹性矩阵,单元应力矩阵
%AREAHKASLODASDISPFP1
%单元面积,总体刚度矩阵,总体荷载向量,结点位移向量,数据文件指针
formatshorte
clear
FP1=fopen('C:
\Users\Administrator\Desktop\input.txt','rt');
NPION=fscanf(FP1,'%d',1);%结点个数(结点编码总数)
NELEM=fscanf(FP1,'%d',1);%单元个数(单元编码总数)
NFORCE=fscanf(FP1,'%d',1);%结点荷载个数
NVFIX=fscanf(FP1,'%d',1);%受约束边界点数
YOUNG=fscanf(FP1,'%e',1);%弹性模量
POISS=fscanf(FP1,'%f',1);%泊松比
THICK=fscanf(FP1,'%d',1);%厚度
LNODS=fscanf(FP1,'%d',[3,NELEM])';%单元定义数组(单元结点号)
COORD=fscanf(FP1,'%f',[2,NPION])';%结点坐标数组
FORCE=fscanf(FP1,'%f',[3,NFORCE])';%结点力数组
FIXED=fscanf(FP1,'%d',[3,inf])';%约束信息数组
%引用所需的全局变量
%globalNPIONNELEMCOORDLNODSYOUNGPOISS
%globalBMATXDMATXSMATXAREA
%生成弹性矩阵D
a=YOUNG/(1-POISS^2);
DMATX(1,1)=1*a;
DMATX(1,2)=POISS*a;
DMATX(2,1)=POISS*a;
DMATX(2,2)=1*a;
DMATX(3,3)=(1-POISS)*a/2;
fori=1:
NELEM;%i为当前所计算的单元号
%计算当前单元的面积
AREA=det([1COORD(LNODS(i,1),1)COORD(LNODS(i,1),2);...
1COORD(LNODS(i,2),1)COORD(LNODS(i,2),2);...
1COORD(LNODS(i,3),1)COORD(LNODS(i,3),2);])/2;
end
%生成应变矩阵B
forj=0:
2
b(j+1)=COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);
c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1);
end
BMATX=[b
(1)0b
(2)0b(3)0;...
0c
(1)0c
(2)0c(3);...
c
(1)b
(1)c
(2)b
(2)c(3)b(3)]/(2*AREA);
SMATX=DMATX*BMATX;%求应力矩阵S=D*B
%所引用的全局变量:
%globalNPIONNELEMNVFIXLNODSASTIFTHICK
%globalBMATXSMATXAREAFIXED
HK=seros(2*NPION,2*NPION);%张成特定大小总刚矩阵并置0
fori=1:
NELEM
EK=BMATX'*SMATX*THICK*AREA;%求解单元刚度矩阵
a=LNODS(i,:
);%临时向量,用来记录当前单元的节点编号
forj=1:
3
fork=1:
3
HK((a(j)*2-1):
a(j)*2,(a(k)*2-1):
a(k)*2)=HK((a(j)*2-1):
a(j)*2,(a(k)*2-1):
a(k)*2)+EK(j*2-1:
j*2,k*2-1:
k*2);
%跟据节点编号对应关系将单元刚度分块叠加到总刚矩阵中
end
end
end
%将约束信息加入总刚(置0置1法)
NUM=1;%计数器(当前已分析的节点数)
i=0;%计数器(当前已处理的约束数)
tmp(NVFIX)=0;%临时存被约束的序号
whilei forj=-1: 0 ifFIXED(NUM,j+3)==1%若发现约束 i=i+1;%计数器自增 tmp(i)=FIXED(NUM)*2+j;%求约束序号 end end NUM=NUM+1; end fori=1: NVFIX HK(1: 2*NPION,tmp(i))=0;%将一约束序号处的总刚列向量清0 HK(tmp(i),1: 2*NPION)=0;%.将一约束序号处的总刚行向量清0 HK(tmp(i),tmp(i))=1;%将行列交叉处的元素置为1 end %所需引用的全局变量 %globalASLODNPIONNFORCEFORCE ASLOD(1: 2*NPION)=0;%张成特定大小的0向量 fori=1: NFORCE ASLOD((FORCE(i,1)*2-1): FORCE(i,1)*2)=FORCE(i,2: 3); end %将有约束处的荷载置为0 NUM=1; i=0; tmp(NVFIX)=0; whilei forj=-1: 0 ifFIXED(NUM,j+3)==1 i=i+1; tmp(i)=FIXED(NUM)*2+j; end end NUM=NUM+1; end fori=1: NVFIX ASLOD(tmp(i))=0; end ASDISP=HK\ASLOD';%计算结点位移向量 %所引用的全局变量 %globalNELEMNPIONSMATXASDISPLNODS ELEDISP(1: 6)=0;%当前单元的结点位移向量 fori=1: NELEM forj=1: 3 ELEDISP(j*2-1: j*2)=ASDISP(LNODS(i,j)*2-1: LNODS(i,j)*2); end STRESS=SMATX*ELEDISP';%求内力 end
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 平面三角形 单元 Matlab 程序