桁架结构MATLAB编程

1、桁架结构% 变量说明% NPOIN NELEM NVFIX NFPOIN NFPRES%总结点数,单元数，约束个数，受力结点数，非结点力数% COORD LNODS YOUNG%结构节点坐标数组，单元定义数组，弹性模量%FPOINFPRESFORCE FIXED%结点力数组，非结点力数组，总体荷载向量，约束信息数组% HKDISP%总体刚度矩阵，结点位移向量%*format short e%设定输出类型clear%清除内存变量FPl=fopen (' 6-6. txt',' rt');%打开初始数据文件%读入控制数据NELEM=fscanf(FPl,'%

2、d',l); %单元数NPOIN=fscanf (FP1,1 %d', 1); %结点数NVFIX=fscanf (FP1,' %d1 1); %约束数NFPOIN=fscanf (FP1,1, 1); %结点荷我数NFPRES=fscanf (FP1,1, 1); %非结点荷载数YOUNG=fscanf(FPl,'%f，l); %弹性模量%读取结构信息LNODS=fscanf(FP1,' %f1, 4, NELEM)'%单元定义：左、右结点号，面积，惯性矩(共计NELEM组)COORD=fscanf (FP1,' %f9, 2, NPO

3、IN)'%坐标：X, y坐标(共计NPOIN组)FPOIN=fscanf(FP1,J , 3, NFPOIN)f%节点力(共计NFPOIN组)：结点号、X方向力(向右正)，% Y方向力(向上正)，FPRES=fscanf (FP1,J , 4, NFPRES)'% 均布力(共计% NFPRES组)：单元号、荷载类型、荷载大小、距离左端长度 FIXED=fscanf(FPl,，%f，, NVFIX)'%约束信息：约束对应的位移编码(共计NVFIX组)HK=zeros (2+NPOIN, 2+NPOIN);%张成总刚矩阵并清零FORCE二zeros (2*NPOIN, 1)

4、;%张成总荷教向量并清零%形成总刚 for i=l:NELEM %对单元个数循环%生成局部单刚(局部坐标)右手坐标系EK=ele.EK(i, LNODS, COORD, YOUNG);T=zbzh(i, LNODS, COORD); % 坐标转换矩阵EKT=T *EK*T;%生成整体单刚(整体坐标系)%组成总刚按2*2子块加入总刚中(共计4块)for j=l:2%对行进行循环按结点号循环N1=LNODS (i, j)*2;% j结点第2个位移的整体编码for k=l:2%对列进行循环一按结点号循环N2=LNODS(i, k)*2; % k结点第2个位移的整体编码HK(N1-1):N1, (N2

5、-1) :N2)=HK(N1-1) :N1, (N2-1):N2) +EKT(j*2-l:j*2, k*2-l:k*2); end end end %由结点力与非结点力生成总荷载向量列阵 for i=l:NFPOIN%对结点荷载个数进行循环Nl=FPOIN(i, 1);%作用荷载的结点号Nl=Nl*2-2;%该结点号对应第一个位移编码- 1for j=l:2FORCE (N1+j) =FORCE (N1+j) +FPOIN (i, j+D ;%取结点荷载 end end%计算由非结点荷载引起的等效结点荷载for i=l：NFPRES %对非结点荷载个数进行循环FO=ele_FPRES (i,

6、FPRES, LNODS, COORD); %计算单元固端力 %对单元局部杆端力要进行坐标转换 T=zbzh(i,LNODS,COORD); % 坐标转换矩阵 FO=T' *F0;ele=FPRES (i, 1); %取荷载所在的单元号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);

7、end%总刚、总荷载进行边界条件处理for j=l：NVFIX %对约束个数进行循环Nl=FIXED(j);HK(l:2*NP0IN,Nl)=0; HK(N1, l:2*NP0IN)=0; HK(N1,N1)=1;%将零位移约束对应的行、列变成零，主元变成1FORCE (Nl)=0;end%-DISP=HKFORCE %方程求解，HK先求逆再与力向量左乘%求结构各个单元内力EDISP=zeros (4, 1); %单元位移列向量清零for i=l:NELEM %对单元个数进行循环for j=l：2%对杆端循环% i单元左右端结点号*2 =该结点的最后一个位移编码 Nl=LNODS(i, j)*

8、2;%取一端的单元位移列向量EDISP(2*j-l:2*j)=DISP(Nl-l:Nl);end%生成局部单刚(局部坐标)右手坐标系EK=ele.EK(i, LNODS, COORD, YOUNG);T=zbzh(i,LNODS, COORD); % 坐标转换矩阵FE=EK*T*EDISP; %计算局部坐标杆端力(由结点位移产生)j=l:NFPRESif FPRES(j, 1) = i %成立时，当前单元上有非结点荷载 FO=ele.FPRES (j, FPRES, LNODS, COORD);%单元固端力FE=FE+F0; %考虑由非结点荷载引起的杆端力end%打印杆端力end FE end

9、%计算桁架单元刚度矩阵函数EK%入口参数：单元号、单元信息数组、结点坐标、弹性模量 %出口参数：局部桁架单元刚度矩阵EK function EK=ele.EK(i> LNODS, COORD, YOUNG)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); % 面积； E=YOUNG;%生成单刚(局部坐标)右手坐标系EK =E*A/L

10、0-E*A/L000000E*A/L00 ；0 ；0 1；return%计算单元固端力函数(正方向：X向右Y向上M逆时针)%入口参数：荷载序号，荷载信息，单元信息，结点坐标%出口参数：单元固端力左右两端的轴力、剪力、弯矩function FO=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=C

11、OORD(NR, 1)-COORD (NL, 1);% x 坐标差dy=COORD (NR, 2) -COORD (NL, 2);% y 坐标差L=sqrt (dx"2+dyA2);%单元长度%计算公式中一些常出现的项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;%单元固端力清零switch FPRES(iFPRES, 2)case 1 %均布荷载FO (2)=-G*C*(2-2*C2+C3)/2. 0;FO (3)=-G*C*C* (6-8*Cl+3*C2)/12. 0;F0(5)=-G*C-F0(2);

12、FO (6) =G*C*C*C1* (4-3*Cl)/12. 0;case 2 %横向集中力FO(2)="G*B1*B2*(L+2*C);FO(3)=-G*C*B1*B1;FO (5)="G*C2* (L+2*D)/L;F0(6)=G*D*C2;case 3 %纵向集中力FO(1)=-G*B1;FO(4)=-G*C1;endreturn%形成第i单元的坐标转换矩阵函数T(4, 4)%入口参数：单元号，单元信息，结点坐标%出口参数：坐标转换矩阵(整体向局部投影)function T=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); %