四边形四结点等参元matlab程序广州大学.docVIP

%---------------四边形四节点等参元 matlab计算程序---------------------------- % 2013年 % 13级建筑与土木工程 PHILDER % B15-404 % 本程序只输出了节点位移、单元中心点的应力 %******************************************************************* clear all; %清楚内存变量 clc; %清屏 format short e %设定数据输出类型 FP1=fopen(input.txt,rt);%打开初始数据文件 NELEM=fscanf(FP1,%d,1); %单元数 NPOIN=fscanf(FP1,%d,1); % 总结点数 NFPOIN=fscanf(FP1,%d,1); %受力结点数 NVFIX=fscanf(FP1,%d,1); %约束结点个数 E=fscanf(FP1,%f,1); %弹性模量 v=fscanf(FP1,%f,1); % 泊松比 t=fscanf(FP1,%f,1); %厚度 LNODS=fscanf(FP1,%f,[4,NELEM]); % 单元定义： 单元结点号（逆时针） COORD=fscanf(FP1,%f,[2,NPOIN]); % 结点号 x,y坐标(整体坐标下) FPOIN=fscanf(FP1,%f,[3,NFPOIN]); % 节点力：结点号、X方向力(向右正)，Y方向力(向上正) FIXED=fscanf(FP1,%d,[3,NVFIX]); %约束信息数组(n,3) n:受约束节点数目, (n,1):约束点号 (n,2),(n,3)约束对应的位移编码 HK=zeros(2*NPOIN ,2*NPOIN); %生成总刚度矩阵并清零 FORCE=zeros(2*NPOIN,1); %生成总荷载列阵并清零 for i=1:NELEM %对单元个数进行循环 EK=ele_EK(i, E, v, LNODS, COORD ,t); %生成单刚矩阵 for j=1:4 %对行循环 N1=LNODS(i,j)*2; %j结点对应的第二个位移编码 for k=1:4 %对列循环 N2=LNODS(i,k)*2; %k结点对应的第二个位移编码 HK(N1-1:N1,N2-1:N2)= HK(N1-1:N1,N2-1:N2)+EK(j*2-1:j*2,k*2-1:k*2); end end end for i=1:NFPOIN %对作用荷载的结点数进行循环 N1=FPOIN(i,1)*2-2; %作用荷载的节点对应的第一个位移编码-1 for j=1:2 %X, Y方向 FORCE(N1+j)=FPOIN(i, j+1); end end for i=1:NVFIX %对约束数进行循环 N1=FIXED(i); %约束对应的位移编码 HK(:,N1)=0;HK(N1,:)=0; HK(N1,N1)=1; FORCE(N1,1)=0; %将零位移约束对应的行、列变成零，主元变成1 end DISP=HK\FORCE; %方程求解，得结点位移 a=E/(1-v*v); D=[a,a*v,0;a*v,a,0;0,0,a*(1-v)/2];%弹性矩阵 B=ele_B(i, -1, -1, COORD,LNODS); %计算B矩阵，并赋值（-1，-1）用以验算 EDISP=zeros(8,1); %生成单元位移列向量，并清零 S=D*B; %计算应力矩阵S for i=1:NELEM %对单元数进行循环 for j=1:4 %对结点数进行循环