planeframe.txt

平面刚架（带弹簧支座）有限元程序

fp1=fopen('data.txt','r');
fp2=fopen('result','w');
%1.读入结构数据，建立累积约束表向量，建立结构刚度矩阵
%1.1 结构参数和弹性模量
[m,c]=fscanf(fp1,'%u',[1]);%杆件数
[nj,c]=fscanf(fp1,'%u',[1]);%节点数
[nrj,c]=fscanf(fp1,'%u',[1]);%约束节点数
[nej,c]=fscanf(fp1,'%u',[1]);%弹性支座节点数
[e,c]=fscanf(fp1,'%e',[1]);%弹性模量

fprintf(fp2,'(1)结构参数和弹性模量\n');
fprintf(fp2,'杆件数  节点数  约束节点数 弹性支座节点数 弹性模量\n');
fprintf(fp2,'%3u  %8u  %8u  %13.3e\n',m,nj,nrj,nej,e);
fprintf(fp2,'\n');
  
%1.2 节点坐标
pc(1:2,1:nj)=0;
for jx=1:nj
[k,c]=fscanf(fp1,'%u',[1]); %节点号
[pc(:,k),c]=fscanf(fp1,'%f',[2]);%x坐标 y坐标
end
fprintf(fp2,'(2)结点坐标\n');
fprintf(fp2,'节点号  x坐标   y坐标\n');
fprintf(fp2,'%3u  %13.3e%13.3e\n',[1:nj;pc(:,1:nj)]);
fprintf(fp2,'\n');


%1.3杆件标号与截面特征
jj(1:m)=0;jk(1:m)=0;ax(1:m)=0;iz(1:m)=0;
for imx=1:m
[k,c]=fscanf(fp1,'%u',[1]);%杆件号
[jj(k),c]=fscanf(fp1,'%u',[1]);%j端点号
[jk(k),c]=fscanf(fp1,'%u',[1]);%k端点号
[ax(k),c]=fscanf(fp1,'%f',[1]);%截面积
[iz(k),c]=fscanf(fp1,'%f\n',[1]);%惯性矩
end

fprintf(fp2,'(3)杆件标号与截面特性\n');
fprintf(fp2,'杆件号  j端点号  k端点号  截面积  截面惯性矩\n');
fprintf(fp2,'%3u  %8u%8u  %13.3e%13.3e\n',[1:m;jj(1:m);jk(1:m);ax(1:m);iz(1:m)]);
fprintf(fp2,'\n');

%1.4节点约束情况
fprintf(fp2,'(4)节点约束情况\n');
fprintf(fp2,'受约束点号  x向约束情况  y向约束情况  z向约束情况\n');
rl(1:3*nj)=0;
for jx=1:nrj
[k,c]=fscanf(fp1,'%u',[1]);%受约束节点号
[rl(3*k-2:3*k),c]=fscanf(fp1,'%f',[3]);%x向约束情况，y向约束情况，z向约束情况
fprintf(fp2,'%5u %12u%12u%12u\n',k,rl(3*k-2:3*k));
end
fprintf(fp2,'\n');
%1.5建立累积约束表向量
crl(1:3*nj)=0;
[n,ne,nr,crl]=CummulatedRestrainedList3(3*nj,rl);

fprintf(fp2,'(5)非约束位移数 弹性约束位移数  约数位移数\n');
fprintf(fp2,'%10u %10u %10u\n',n,ne,nr);
fprintf(fp2,'\n');
%1.55读入弹性支座刚度系数
sr(1:3*nj,1:1)=0;
fprintf(fp2,'(5.5)弹性支座刚度系数\n');
fprintf(fp2,'节点号 x向线位移  y向线位移 z向角线位移\n');
for jx=1:nej
[kr,c]=fscanf(fp1,'%u',[1]);%发生约束位移的节点号
[sr(3*kr-2:3*kr),c]=fscanf(fp1,'%e',[3]);%x方向线位移，y方向线位移，z向角线位移
fprintf(fp2,'%3u%14.3e%14.3e%14.3e\n',kr,sr(3*kr-2:3*kr));
end
fprintf(fp2,'\n');
%对应非约束位移，约束位移，将支座位移向量sr重新排列
sr(crl(1:3*nj))=sr(1:3*nj);

%1.6装配总结点刚度矩阵
sj(1:3*nj,1:3*nj)=0;
for im=1:m
smd=PlaneFrameMemberGlobalStiffnessMatrix(pc(:,jj(im)),pc(:,jk(im)),ax(im),iz(im),e);
j=[3*jj(im)-2:3*jj(im)   3*jk(im)-2:3*jk(im)];%计算杆端位移对应的总节点位移编号
sj(j(1:6),j(1:6))=sj(j(1:6),j(1:6))+smd(1:6,1:6);
end
%1.7对应非约束位移，约束位移，将总节点刚度矩阵sj重新排列
s(crl(1:3*nj),crl(1:3*nj))=sj(1:3*nj,1:3*nj);


%2读入荷载数据，构建约数杆端力，构建综合节点荷载
%2.1读入荷载节点数，集中荷载杆件数，分布荷载杆件数
[nlj,c]=fscanf(fp1,'%u',[1]);%荷载节点数
[nla,c]=fscanf(fp1,'%u',[1]);%集中荷载杆件数
[nlq,c]=fscanf(fp1,'%u',[1]);%分布荷载杆件数

fprintf(fp2,'(6)荷载节点数  集中荷载杆件数  分布荷载杆件数\n');
fprintf(fp2,'%10u  %10u  %10u\n',nlj,nla,nlq);
fprintf(fp2,'\n');

%2.2结点荷载
fprintf(fp2,'(7)结点荷载\n');
fprintf(fp2,'节点号  x向线力  y向线力  z向力偶\n');
ac(1:3*nj,1:1)=0;
for jx=1:nlj
[k,c]=fscanf(fp1,'%u',[1]);%受荷载节点号
[aj,c]=fscanf(fp1,'%f',[3]);%x向线力 y向线力 z向力偶
fprintf(fp2,'%3u  %13.3e%13.3e%13.3e\n',k,aj);
j=[3*k-2:3*k];
ac(j(1:3))=ac(j(1:3))+aj(1:3);
end
fprintf(fp2,'\n');


%2.3杆件上集中荷载
fprintf(fp2,'(8)杆件上集中荷载\n');
fprintf(fp2,'杆件号  x向线力  y向线力  z向力偶  作用点距杆j端距离\n');
aml(1:6,1:m)=0;
for imx=1:nla
[ia,c]=fscanf(fp1,'%u',[1]);%杆件号
[as,c]=fscanf(fp1,'%f',[3]);%x向线力 y向线力 z向力偶
[xa,c]=fscanf(fp1,'%f',[1]);%集中荷载杆端距
fprintf(fp2,'%3u  %13.3e%13.3e%13.3e%13.3e\n',ia,as,xa);
amx=PlaneFrameMemberConcentratedLoad(pc(:,jj(ia)),pc(:,jk(ia)),as,xa);
aml(:,ia)=aml(:,ia)+amx(:);%累加约束杆端力
end
fprintf(fp2,'\n');


%2.4 杆件上分布荷载
fprintf(fp2,'(9)杆件上分布荷载\n');
fprintf(fp2,'杆件号 起点x向集度  终点x向集度  起点y向集度  终点y向集度  起点距j端距离  终点距j端距离\n');
for imx=1:nlq
[iq,c]=fscanf(fp1,'%u',[1]);%杆件号
[q,c]=fscanf(fp1,'%f',[4]);%起点x向集度 终点x向集度 起点y向集度 终点y向集度
[xq,c]=fscanf(fp1,'%f',[2]);%起点距杆j端距离，终点距杆j端距离
fprintf(fp2,'%3u %14.3e%14.3e%14.3e%14.3e%14.3e%14.3e\n',iq,q,xq);
amx=PlaneFrameMemberDistributionLoad(pc(:,jj(iq)),pc(:,jk(iq)),q,xq);
aml(:,iq)=aml(:,iq)+amx(:);%累加约束杆端力
end
fprintf(fp2,'\n');
fclose(fp1);%关闭输入数据文件

fprintf(fp2,'(10)约束杆端力\n');
fprintf(fp2,'杆件号  j/k端号  xm向线约束力  ym向线约束力  zm向约束弯矩\n');
fprintf(fp2,'%3u %6u %14.3e%14.3e%14.3e\n    %6u %14.3e%14.3e%14.3e\n',...
[1:m;jj(1:m);aml(1:3,1:m);jk(1:m);aml(4:6,1:m)]);
fprintf(fp2,'\n');


%2.5计算等效节点荷载 
ae(1:3*nj,1:1)=0;
for im=1:m
RT=PlaneFrameRotationMatrix(pc(:,jj(im)),pc(:,jk(im)));
j=[3*jj(im)-2:3*jj(im) 3*jk(im)-2:3*jk(im)];
ae(j(1:6))=ae(j(1:6))-RT'*aml(1:6,im);
end
fprintf(fp2,'(11)等效节点荷载\n');
fprintf(fp2,'节点号  x向线力  y向线力  z向力偶\n');
for ja=1:nj
fprintf(fp2,'%3u%13.3e%13.3e%13.3e\n',ja,ae(3*ja-2:3*ja));
end
fprintf(fp2,'\n');
%2.6计算综合节点荷载
ac(1:3*nj,1:1)=ac(1:3*nj,1:1)+ae(1:3*nj,1:1);
fprintf(fp2,'(12)综合节点荷载\n');
fprintf(fp2,'节点号  x向线力  y向线力  z向力偶\n');
for ja=1:nj
  fprintf(fp2,'%3u %13.3e%13.3e%13.3e\n',ja,ac(3*ja-2:3*ja));
end
fprintf(fp2,'\n');
%2.7对应非约束位移,约束位移，并将综合节点荷载ac重新排列
ac(crl(1:3*nj))=ac(1:3*nj);
%3计算非约束位移，支座反力，并构建总节点位移向量
[d,r]=StiffnessMatrixEquationSolutionWithElasticRestraint(n,ne,nr,s,ac,sr,crl);
fprintf(fp2,'(13)节点位移与支座反力\n');
fprintf(fp2,'节点号   x向线位移  y向线位移  z向角位移  x向支座反力  y向支座反力  z向支座反力矩\n');
for je=1:nj
  fprintf(fp2,'%3u %13.3e%13.3e%13.3e%13.3e%13.3e%13.3e\n',je,d(3*je-2:3*je),r(3*je-2:3*je));
end
fprintf(fp2,'\n'); 
 %4计算最终杆端力
for im=1:m
RT=PlaneFrameRotationMatrix(pc(:,jj(im)),pc(:,jk(im)));
sm=PlaneFrameMemberLocalStiffnessMatrix(pc(:,jj(im)),pc(:,jk(im)),ax(im),iz(im),e);
j=[3*jj(im)-2:3*jj(im),3*jk(im)-2:3*jk(im)];
aml(:,im)=aml(:,im)+sm*RT*d(j(1:6));
end
fprintf(fp2,'(14)最终杆端力\n');
fprintf(fp2,'杆件号j/k端号   xm向轴力   ym向剪力    zm向弯矩\n');
fprintf(fp2,'%3u  %7u  %14.3e%14.3e%14.3e\n     %7u  %14.3e%14.3e%14.3e\n',...
[1:m;jj(1:m);aml(1:3,1:m);jk(1:m);aml(4:6,1:m)]);
fclose(fp2);