平面三角形单元Matlab程序

%变量说明

%NPOIN NELEM NVFIX NFORCE NNODE

%总结点数，单元数，受约束边界点数，结点力数, 单元结点数

%COORD LNODS YOUNG POISS THICK %结构结点坐标数组，单元定义数组，弹性模量，泊松比，厚度

%FORCE FIXED BMA TX 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 BMA TX DMATX SMA TX 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

BMA TX=[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);

SMA TX=DMATX*BMATX; %求应力矩阵S=D*B

% 所引用的全局变量：

%global NPION NELEM NVFIX LNODS ASTIF THICK

%global BMA TX SMA TX AREA FIXED

HK=seros(2*NPION,2*NPION); %张成特定大小总刚矩阵并置0

for i=1:NELEM

EK=BMA TX'*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 SMA TX ASDISP LNODS

ELEDISP(1:6)=0; %当前单元的结点位移向量for i=1:NELEM

for j=1:3