2014基本平面刚架MATLAB程序2014

%平面刚架MATLAB程序%2003.9.162007.2.282008.4.12009.102011.102013.92014.09%*************************************************%变量说明%NPOINNELEMNVFIXNFPOINNFPRES%总结点数,单元数,约束个数,受力结点数,非结点力数%COORDLNODSYOUNG%结构节点坐标数组,单元定义数组,弹性模量%FPOINFPRESFORCEFIXED%结点力数组，非结点力数组，总体荷载向量,约束信息数组%HKDISP%总体刚度矩阵，结点位移向量%**************************************************formatshorte%设定输出类型clear%清除内存变量FP1=fopen('6-6.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%-------------------------------------------------------ele_FPRES.m%计算单元固端力函数（正方向：X向右Y向上M逆时针）%入口参数：荷载序号，荷载信息，单元信息，结点坐标%出口参数：单元固端力——左右两端的轴力、剪力、弯矩functionF0=ele_FPRES(iFPRES,FPRES,LNODS,COORD)ele=FPRES(iFPRES,1);%取荷载所在的单元号G=FPRES(iFPRES,3);%单元荷载大小C=FPRES(iFPRES,4);%单元荷载与左端距离NL=LNODS(ele,1);NR=LNODS(ele,2);%单元的左右结点号dx=COORD(NR,1)-COORD(NL,1);%x坐标差dy=COORD(NR,2)-COORD(NL,2);%y坐标差L=sqrt(dx^2+dy^2);%单元长度%计算公式中一些常出现的项D=L-C;C1=C/L;C2=C1*C1;C3=C1*C2;B1=D/L;B2=B1/L;F0=[0;0;0;0;0;0];%单元固端力清零switchFPRES(iFPRES,2)case1%均布荷载F0(2)=-G*C*(2-2*C2+C3)/2.0;F0(3)=-G*C*C*(6-8*C1+3*C2)/12.0;F0(5)=-G*C-F0(2);F0(6)=G*C*C*C1*(4-3*C1)/12.0;case2%横向集中力F0(2)=-G*B1*B2*(L+2*C);F0(3)=-G*C*B1*B1;F0(5)=-G*C2*(L+2*D)/L;F0(6)=G*D*C2;case3%纵向集中力F0(1)=-G*B1;F0(4)=-G*C1;endreturnele_EK.m%计算单元刚度矩阵函数EK%入口参数：单元号、单元信息数组、结点坐标、弹性模量%出口参数：局部单元刚度矩阵EKfunctionEK=ele_EK(i,LNODS,COORD,E)NL=LNODS(i,1);NR=LNODS(i,2);%左右结点号dx=COORD(NR,1)-COORD(NL,1);%x坐标差dy=COORD(NR,2)-COORD(NL,2);%y坐标差L=sqrt(dx^2+dy^2);%单元长度A=LNODS(i,3);I=LNODS(i,4);%面积；惯性矩%生成单刚(局部坐标)右手坐标系EK=[E*A/L00-E*A/L00;...012*E*I/L^36*E*I/L^20-12*E*I/L^36*E*I/L^2;...06*E*I/L^24*E*I/L0-6*E*I/L^22*E*I/L;...-E*A/L00E*A/L00;...0-12*E*I/L^3-6*E*I/L^2012*E*I/L^3-6*E*I/L^2;...06*E*I/L^22*E*I/L0-6*E*I/L^24*E*I/L];return%--------------------------------------------------------zbzh.m%形成第i单元的坐标转换矩阵函数T(6,6)%入口参数：单元号，单元信息，结点坐标%出口参数：坐标转换矩阵（整体向局部投影）functionT=zbzh(i,LNODS,COORD)NL=LNODS(i,1);%左结点号NR=LNODS(i,2);%右结点号dx=COORD(NR,1)-COORD(NL,1);%x坐标差dy=COORD(NR,2)-COORD(NL,2);%y坐标差L=sqrt(dx^2+dy^2);%单元长度c=dx/L;%cosa(与x轴夹角余弦）s=dy/L;%sinaT=[cs0000;...-sc0000;...001000;...000cs0;...000-sc0;...000001];return