平面三角形单元Matlab程序
%变量说明
%NPOIN NELEM NVFIX NFORCE NNODE
%总结点数，单元数，受约束边界点数，结点力数, 单元结点数
%COORD LNODS YOUNG POISS THICK %结构结点坐标数组，单元定义数组，弹性模量，泊松比，厚度
%FORCE FIXED BMATX DMATX SMATX %节点力数组，约束信息数组，单元应变矩阵，单元弹性矩阵，单元应力矩阵%AREA HK ASLOD ASDISP FP1
%单元面积，总体刚度矩阵，总体荷载向量，结点位移向量，数据文件指针format short e
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])'; %约束信息数组
%引用所需的全局变量
%global NPION NELEM COORD LNODS YOUNG POISS %global BMATX DMATX SMATX AREA
%生成弹性矩阵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;
for i=1:NELEM; %i为当前所计算的单元号
%计算当前单元的面积
AREA=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);...
1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);...
1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2);])/2; end
%生成应变矩阵B
for j=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) 0 b(2) 0 b(3) 0;...
0 c(1) 0 c(2) 0 c(3);...
c(1) b(1) c(2) b(2) c(3) b(3)]/(2*AREA);
SMATX=DMATX*BMATX; %求应力矩阵S=D*B
% 所引用的全局变量：
%global NPION NELEM NVFIX LNODS ASTIF THICK
%global BMATX SMATX AREA FIXED
HK=seros(2*NPION,2*NPION); %张成特定大小总刚矩阵并置0
for i=1:NELEM
EK=BMATX'*SMATX*THICK*AREA; %求解单元刚度矩阵
a=LNODS(i,:); %临时向量,用来记录当前单元的节点编号
for j=1:3
for k=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; %临时存被约束的序号
while i for j=-1:0 if FIXED(NUM,j+3)==1 %若发现约束 i=i+1; %计数器自增 tmp(i)=FIXED(NUM)*2+j; %求约束序号 end end NUM=NUM+1; end for i=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 %所需引用的全局变量 %global ASLOD NPION NFORCE FORCE ASLOD(1:2*NPION)=0; %张成特定大小的0向量for i=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; while i for j=-1:0 if FIXED(NUM,j+3)==1 i=i+1; tmp(i)=FIXED(NUM)*2+j; end end NUM=NUM+1; end for i=1:NVFIX ASLOD(tmp(i))=0; end ASDISP=HK\ASLOD'; %计算结点位移向量 %所引用的全局变量 %global NELEM NPION SMATX ASDISP LNODS ELEDISP(1:6)=0; %当前单元的结点位移向量for i=1:NELEM for j=1:3 ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2); end STRESS=SMATX*ELEDISP'; %求内力end THANKS !!! 致力为企业和个人提供合同协议，策划案计划书，学习课件等等 打造全网一站式需求 欢迎您的下载，资料仅供参考