2. Matlab语言程序代码
%******************************************************************************
%初始化及数据调用
format short e %设定输出类型 clear %清除内存变量
FP1=fopen('','rt'); %打开输入数据文件，读入控制数据
NELEM=fscanf(FP1,'%d',1), %单元个数（单元编码总数） NPION=fscanf(FP1,'%d',1), %结点个数（结点编码总数） NVFIX=fscanf(FP1,'%d',1) %受约束边界点数 NFORCE=fscanf(FP1,'%d',1), %结点荷载个数 YOUNG=fscanf(FP1,'%e',1), %弹性模量 POISS=fscanf(FP1,'%f',1), %泊松比 THICK=fscanf(FP1,'%f',1) %厚度
LNODS=fscanf(FP1,'%d',[3,NELEM])' %单元定义数组（单元结点号）
%相应为单元结点号（编码）、按逆时针顺序输入
COORD=fscanf(FP1,'%f',[2,NPION])' %结点坐标数组
%坐标：x，y坐标（共NPOIN组）
FORCE=fscanf(FP1,'%f',[3,NFORCE])' %结点力数组（受力结点编号, x方向,y方向） FIXED=fscanf(FP1,'%d',[3,NVFIX])' %约束信息（约束点，x约束，y约束） %有约束为1，无约束为0
%*****************************************************************************
%生成单元刚度矩阵并组成总体刚度矩阵
ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0
%***************************************************************************** for i=1:NELEM
%生成弹性矩阵D
D= [1 POISS 0; POISS 1 0;
0 0 (1-POISS)/2]*YOUNG/(1-POISS^2)
%*****************************************************************************
%计算当前单元的面积A
A=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
%*****************************************************************************
%生成应变矩阵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
B=[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*A); B1( :,:,i)=B;
%*****************************************************************************
%求应力矩阵S=D*B
S=D*B;
ESTIF=B'*S*THICK*A; %求解单元刚度矩阵
a=LNODS(i,:); %临时向量,用来记录当前单元的节点编号 for j=1:3 for k=1:3
ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)…
=ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2); %根据节点编号对应关系将单元刚度分块叠加到总刚
%度矩阵中
end end end
%*****************************************************************************
%将约束信息加入总体刚度矩阵（对角元素改一法）
for i=1:NVFIX
if FIXED(i,2)==1
ASTIF(:,(FIXED(i,1)*2-1))=0; %一列为零 ASTIF((FIXED(i,1)*2-1),:)=0; %一行为零 ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1; %对角元素为1 end
if FIXED(i,3)==1
ASTIF( :,FIXED(i,1)*2)=0; %一列为零 ASTIF(FIXED(i,1)*2,:)=0; %一行为零 ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; %对角元素为1 end end
%*****************************************************************************
%生成荷载向量
ASLOD(1:2*NPION)=0; %总体荷载向量置零 for i=1:NFORCE
ASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3); end
%*****************************************************************************
%求解内力
ASDISP=ASTIF\\ASLOD' %计算节点位移向量 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 i
STRESS=D*B1(:, :, i)*ELEDISP' %求内力 end
fclose(FP1); %关闭数据文件
五、 程序运行结果
NELEM = 4 NPION = 6 NVFIX = 2 NFORCE = 1 YOUNG = 0 POISS =
THICK =
LNODS =
1 2 6 2 3 4 2 4 5 2 5 6 COORD =
0 0 1 0 2 0 2 1 1 1 0 1 FORCE =
4 0 -150 FIXED =
1 1 1 6 1 1 D =
+008 +007 0 +007 +008 0
0 0 +007 A =
D =
+008 +007 0 +007 +008 0
0 0 +007 A = D =
+008 +007 0 +007 +008 0
0 0 +007 A = D =
+008 +007 0 +007 +008 0
0 0 +007 A =
ASDISP =
0 0
0 0 （说明：以上12个ASDISP输出结果数据表示节点位移向量，即各节点分别在X，Y方向上产生的位移。） i = 1 STRESS = +005 +004 +005 i = 2 STRESS = +004
相关推荐：
loading