基本平面刚架MATLAB程序

%平面刚架MATLAB程序%2003.9.162007.2.282008.4.12009.102011.102013.92014.09%*************************************************%变量说明%NPOINNELEMNVFIXNFPOINNFPRES%总结点数,单元数,约束个数,受力结点数,非结点力数%COORDLNODSYOUNG%结构节点坐标数组,单元定义数组,弹性模量A%FPOINFPRESFORCEFIXED%结点力数组，非结点力数组，总体荷载向量,约束信息数组%HKDISP%总体刚度矩阵，结点位移向量%**************************************************formatshorte%设定输出类型clear%清除内存变量FP1=fopen('XX.txt','rt')%打开初始数据文件%读入控制数据NELEM=fscanf(FP1,'%d',1);%单元数NPOIN=fscanf(FP1,'%d',1);%结点数NVFIX=fscanf(FP1,'%d',1);%约束数NFPOIN=fscanf(FP1,'%d',1);%作用荷载的结点个数NFPRES=fscanf(FP1,'%d',1);%非结点荷载数YOUNG=fscanf(FP1,'%f',1);%弹性模量%读取结构信息LNODS=fscanf(FP1,'%f',[4,NELEM])'%单元定义：左、右结点号，面积，惯性矩（共计NELEM组）COORD=fscanf(FP1,'%f',[2,NPOIN])'%坐标：x,y坐标（共计NPOIN组）FPOIN=fscanf(FP1,'%f',[4,NFPOIN])'%节点力（共计NFPOIN组）：受力结点号、X方向力（向右正），%Y方向力（向上正），M力偶（逆时针正）FPRES=fscanf(FP1,'%f',[4,NFPRES])'%均布力（共计%NFPRES组）：单元号、荷载类型、荷载大小、距离左端长度%荷载类型1-均布荷载2-横向集中力3-纵向集中力FIXED=fscanf(FP1,'%f',NVFIX)'%约束信息：约束对应的位移编码（共计NVFIX组）%---------------------------------------------------------HK=zeros(3*NPOIN,3*NPOIN);%张成总刚矩阵并清零FORCE=zeros(3*NPOIN,1);%张成总荷载向量并清零%形成总刚fori=1:NELEM%对单元个数循环%生成局部单刚(局部坐标)右手坐标系EK=ele_EK(i,LNODS,COORD,YOUNG);T=zbzh(i,LNODS,COORD);%坐标转换矩阵EKT=T'*EK*T；%生成整体单刚(整体坐标系)%组成总刚按3*3子块加入总刚中（共计4块)forj=1:2%对行进行循环---按结点号循环N1=LNODS(i,j)*3;%j结点第3个位移的整体编码fork=1:2%对列进行循环---按结点号循环N2=LNODS(i,k)*3;%k结点第3个位移的整体编码HK((N1-2):N1,(N2-2):N2)=HK((N1-2):N1,(N2-2):N2)...+EKT(j*3-2:j*3,k*3-2:k*3);End%单刚3x3子块叠加到总刚中endend%由结点力与非结点力生成总荷载向量列阵fori=1:NFPOIN%对结点荷载个数进行循环N1=FPOIN(i,1);%作用荷载的结点号N1=N1*3-3;%该结点号对应第一个位移编码-1forj=1:3FORCE(N1+j)=FORCE(N1+j)+FPOIN(i,j+1);%取结点荷载endend%计算由非结点荷载引起的等效结点荷载fori=1:NFPRES%对非结点荷载个数进行循环F0=ele_FPRES(i,FPRES,LNODS,COORD);%计算单元固端力%对单元局部杆端力要进行坐标转换ele=FPRES(i,1);%取荷载所在的单元号T=zbzh(ele,LNODS,COORD);%坐标转换矩阵F0=T’*F0;NL=LNODS(ele,1);NR=LNODS(ele,2);%单元的左右结点号%将单元固端力变成等效结点荷载（注意固端力与等效结点荷载符号相反）FORCE((3*NL-2):3*NL)=FORCE((3*NL-2):3*NL)-F0(1:3);FORCE((3*NR-2):3*NR)=FORCE((3*NR-2):3*NR)-F0(4:6);end%总刚、总荷载进行边界条件处理forj=1:NVFIX%对约束个数进行循环N1=FIXED(j);HK(1:3*NPOIN,N1)=0;HK(N1,1:3*NPOIN)=0;HK(N1,N1)=1;%将零位移约束对应的行、列变成零，主元变成1FORCE(N1)=0;end%---------------------------------------------------------DISP=HK\FORCE%方程求解，HK先求逆再与力向量左乘%求结构各个单元内力EDISP=zeros(6,1);%单元位移列向量清零fori=1:NELEM%对单元个数进行循环forj=1:2%对杆端循环%i单元左右端结点号*3=该结点的最后一个位移编码N1=LNODS(i,j)*3;%取一端的单元位移列向量EDISP(3*j-2:3*j)=DISP(N1-2:N1);end%生成局部单刚(局部坐标)右手坐标系EK=ele_EK(i,LNODS,COORD,YOUNG);T=zbzh(i,LNODS,COORD);%坐标转换矩阵FE=EK*T*EDISP;%计算局部坐标杆端力（由结点位移产生）forj=1:NFPRESifFPRES(j,1)==i%成立时，当前单元上有非结点荷载F0=ele_FPRES(j,FPRES,LNODS,COORD);%单元固端力FE=FE+F0;%考虑由非结点荷载引起的杆端力endendFE%打印杆端力end