资源描述:
《基本平面刚架MATLAB程序2011.doc》由会员上传分享,免费在线阅读,更多相关内容在教育资源-天天文库。
1、%平面刚架MATLAB程序%2003.9.162007.2.282008.4.12009.102011.10%*************************************************%变量说明%NPOINNELEMNVFIXNFPOINNFPRES%总结点数,单元数,约束个数,受力结点数,非结点力数%COORDLNODSYOUNG%结构节点坐标数组,单元定义数组,弹性模量%FPOINFPRESFORCEFIXED%结点力数组,非结点力数组,总体荷载向量,约束信息数组%HKDISP%总体刚度矩阵,结点位移向量
2、%**************************************************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',
3、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
4、,NFPRES])'%均布力(共计%NFPRES组):单元号、荷载类型、荷载大小、距离左端长度FIXED=fscanf(FP1,'%f',NVFIX)'%约束信息:约束对应的位移编码(共计NVFIX组)%---------------------------------------------------------HK=zeros(3*NPOIN,3*NPOIN);%张成总刚矩阵并清零FORCE=zeros(3*NPOIN,1);%张成总荷载向量并清零%形成总刚fori=1:NELEM%对单元个数循环%生成局部单刚(局部坐标)右手
5、坐标系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)...
6、+EKT(j*3-2:j*3,k*3-2:k*3);endendend%由结点力与非结点力生成总荷载向量列阵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);%计算单元固
7、端力%对单元局部杆端力要进行坐标转换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%总刚、总荷载进行
8、边界条件处理forj=1:NVFIX%对约束个数进行循环N1=FIXED(j);HK(1:3*NPOIN,N1)=0;HK(N1,1:3*NPOIN)=0;HK(N1,N1)=1;%将零位移约束对应的行、列变成零,主元变成1FORCE