matlab.zip_dooruxc_matlab 刚架_结构位移

clear; %清除内存变量 ss1=input('give a data filename:','s'); fp1=fopen(ss1,'r'); [msg,errn]=ferror(fp1); ss2=input('give a result filename:','s'); fp2=fopen(ss2,'w'); [msg,errn]=ferror(fp2); %1.读入结构数据、建立累积约束表向量、建立结构刚度矩阵 %1.1.结构参数和弹性模量 [m,c]=fscanf(fp1,'%u',[1]); %杆件数 [nj,c]=fscanf(fp1,'%u',[1]); %结点数 [nrj,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,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]); 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,nr,crl]=CumulatedRestrainedList(3*nj,rl); fprintf(fp2,'(5)非约束位移数 约束位移数\n'); fprintf(fp2,'%10u %10u\n',n,nr); fprintf(fp2,'\n'); %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); ac(3*k-2:3*k)=ac(3*k-2:3*k)+aj(1:3); end fprintf(fp2,'\n'); %杆件上集中荷载 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.计算等效节点荷载，累加到综合节点荷载中 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)] ac(j(1:6))=ac(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,ac(3*ja-2:3*ja)); end fprintf(fp2,'\n'); %2.6.对应非约束位移，约束位移，将综合节点荷载ac重新排列 ac(crl(1:3*nj))=ac(1:3*nj); %3.计算非约束位移、支座反力，并构建总节点位移向量 [d,r]=StiffnessMatrixEquationSolution(n,nr,s,ac,crl); fprintf(fp2,'(12)节点位移与支座反力\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'); %计算最终杆端力 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,'(13)最终杆端力\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); %关闭计算结果数据文件