2024基本平面刚架各种荷载MATLAB程序.docx

《2024基本平面刚架各种荷载MATLAB程序.docx》由会员分享，可在线阅读，更多相关《2024基本平面刚架各种荷载MATLAB程序.docx（5页珍藏版）》请在优知文库上搜索。

1、%平面刚架MATLAB程序%2024.9.162024.4.12024.102024.102024.92024.092024.03%变说明%NPOINNELEMNVFIXNFPOINNFPRES%总结点数,单元数，约束个数，受力结点数，非结点力数%COORDLNODSYOUNG%结构节点坐标数组，单元定义数组，弹性模%FPOINFPRESFORCEFIXED%结点力数组，非结点力数组，总体荷载向,约束信息数组%HKDISP%总体刚度矩阵，结点位移向量formatshorte%设定输出类型clear%清除内存变FPl=fopenC6-6.txf,rt）%打开初始数据文件%读入限制数据%单元数%结

2、点数%约束数%作用荷载的结点个数%非结点荷载数%弹性模NELEM=fscanf(FPl,%d,l);NPOIN=fscanf(FPl,%d,l);NVFIX=fccanf(FPl；%d,l);NFPOIN=fscanf(FPl,%dl);NFPRES=fscanf(FPl,%dl);YOUNG=fscanf(FPl,%,l);%读取结构信息1.NODS=fscanf（FP1,%,6,NELEM）%单元定义：左、右结点号，面积，惯性矩，线膨胀系数，截面高度（共计NELEM组）COORD=fscanf（FPl,%,2,NPOIND,%坐标：X,y坐标（共计NPOIN组）FPOIN=fccanf（F

3、Pl,%f,4,NFPOIN）,%节点力（共计NFPOIN组受力结点号、X方向力（向右正%Y方向力（向上正），M力偶（逆时针正）FPRES=fscanf（FPl,%f,7,NFPRES）*%均布力（共计%NFPRES组单元号、荷或类型、荷载大小、距离左端长度，温差二（下端-上端）梯形上边。下边（改）%荷载类型1均布荷载2.横向集中力3-纵向集中力4三角形荷载5.温度荷裁6.梯形荷裁FIXED=fscanf（FPl,%,NVFIX）,%约束信息：约束对应的位移编码（共计NVFlX组）%HK=zeros（3*NPOIN,3*NPOIN）;%张成总刚矩阵并清零FORCE=zeros（3*NPOIN,

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

5、(N2-2):N2).+EKT(j*3-2:j*3,k*3-2：k*3);%单刚3x3子块叠加到总刚中endend%由结点力与非结点力生成总荷栽向列阵fori=l:NFPOIN%对结点荷教个数进行循环Nl=FPOIN(i,l);%作用荷裁的结点号NI=NI*3-3;%该结点号对应第一个位移编码-1forj=l:3FORCE(Nl+j)=FORCE(Nl+j)+FPOIN(i9j+l%取结点荷教endend%计算由非结点荷载引起的等效结点荷载fori=hNFPRES%对非结点荷载个数进行循环FO=ele.FPRES(i,FPRES,LNODS,COORD,YOUNG);%计算单元固端力%对多元局

6、部杆端力要进行坐标转换ele=FPRES(i,l);%取荷载所在的单元号T=zbzh(ele,LNODS,COORD);%坐标转换矩阵FO=T,*FO;NL=LNODS(ele,l);NR=LNODS(ele,2);%单元的左右结点号%将单元固端力变成等效结点荷载(留意固端力与等效结点荷载符号相反)FORCE(3*NL-2):3*NL)=FORCE(3*NL-2):3*NL)-FO(1:3);FORCE(3*NR-2):3*NR)=FORCE(3*NR-2):3*NR)-FO(4:6);end%总刚、总荷载进行边界条件处理forj=lNVFIX%对约束个数进行循环Nl=FIXED(j);HK(

7、1:3*NPOIN,N1)=0;HK(NlJ:3*NPOIN)=0;HK(Nl5Nl)=I;%将零位移约束对应的行、列变成零，主元变成1FORCE(Nl)=O;end%DISP=HKFORCE%方程求解，HK先求逆再与力向左乘%求结构各个单元内力EDISP=zeros(6,l);%单元位移列向清零fori=hNELEM%对单元个数进行循环forj=l:2%对杆端循环%i单元左右端结点号*3=该结点的最终一个位移编码Nl=LNODS(iJ)*3;%取一端的单元位移列向EDISP(3*j-2:3*j)=DISP(Nl-2:Nl);end%生成局部单刚(局部坐标)右手坐标系EK=ele_EK(i,L

8、NODS,COORD,YOUNG);T=zbzh(i,LNODS,COORD);%坐标转换矩阵FE=EK*T*EDISP;%计算局部坐标杆端力(由结点位移产生)forj=l:NFPRESifFPRES(j,l)=i%成立时，当前单元上有非结点荷载FO=ele_FPRES(J,FPRES,LNODS,COORD,YOUNG);%单元固端力FE=FE+FO;%考虑由非结点荷裁引起的杆端力endendFE%打印杆端力endend%e!e.FPRES.m%计算单元固端力函数(正方向：X向右Y向上M逆时针)%人口参数：荷载序号，荷载信息，单元信息，结点坐标%出口参数：单元固端力左右两端的轴力、剪力、鸾矩

9、functionFO=ele_FPRES(iFPRES,FPRES,LNODS,COORD,E)ele=FPRES(iFPRES,l);%取荷载所在的单元号G=FPRES(iFPRES,3);%单元荷裁大小C=FPRESCiFPRES,4);%单元荷鼓与左端距离W=FPRES(iFPRES,5);%单元下上端温差S=FPRES(iFPRES,6);%梯形长边荷载X=FPRESCiFPRES,6);%梯形短边荷裁H=LNODS(eIe,6);%单元截面高度P=LNODS(ele,5);%单元裁面线膨胀系数NL=LNODS(ele,l);NR=LNODS(ele,2);%单元的左右结点号dx=CO

10、ORD(NR,l)-COORD(NL,l);%x坐标差dy=COORD(NR,2)-COORD(NL,2);%y坐标差1.=Sqrt(dx2+dy人2);%单元长度%计算公式中一座常出现的项D=L-C;C1=CL;C2=C1*C1;C3=C1*C2;Bl=DZL;B2=B1L;F0=0;0;();();0;();%单元固端力清零switchFPRES(iFPRES,2)case 1 %均布荷载F0(2)=-G*C*(2-2*C2+C3)20;F0(3)=-G*C*C*(6-8*Cl+3*C2)12.0;F0(5)=-G*C-F0(2);FO(6)=G*C*C*C1*(4-3*C1)12.O;c

11、ase 2 %横向集中力FO(2)=-G*B1*B2*(L+2*C);Fo=-G*C*Bl*Bl;F0(5)=-G*C2*(L+2*D)L;F0(6)=G*D*C2;case 3 %纵向集中力FO(1)=-G*B1;FO(4)=-G*C1;case 4 %三角形荷载F0(2)=-7*G*L20;F0(3)=G*L220;F0(5)=3*G*L20;F0(6)=G*L*L30;case 5 %温度荷载F0(3)=E*I*W*PH;F0(6)=-E*I*W*PH;case 6 %梯形荷载(改)Fl(2)=-X*C*(2-2*C2+C3)20;Fl(3)=-X*C*C*(6-8*Cl+3*C2)12

12、0;F2(5)=-X*C-F0(2);F2(6)=X*C*C*C1*(4-3*C1)12O;F3(2)=-7*(S-G)*L20;F3(3)=(S-G)*La220;F4(5)=3*(S-G)*L20;F4(6)=(S-G)*L*L30;F0(2)=Fl(2)+F3(2);F0(3)=Fl(3)+F3(3);F0(5)=F2(5)+F4(5);F0(6)=F2(6)+F4(6);endreturnele_EK.m%计算单元刚度矩阵函数EK%人口参数：单元号、单元信息数组、结点坐标、弹性模量%出口参数：局部单元刚度矩阵EKfunctionEK=ele.EK(i,LNODS,COORD,E)NL=

13、LNODS(i,l);NR=LNODS(i,2);%左右结点号dx=COORD(NR,l)-COORD(NL,l);%x坐标差EK=E*AL00-E*AL0012*E*IL36*E*IL20-12*E*IL306*E*IL24*E*IL0-6*E*II2-E*AL00E*AL00-12*E*IL3-6*E*IL2012*E*IL30;.6*E*IL2;,2*E*IL;.0;.-6*E*IL2;.dy=COORD(NR,2)-COORD(NL,2);%y坐标差1.=sqrt(dxA2+dyA2);%单元长度A=LNODS(i,3);I=LNODS(i,4);%面积；惯性矩生成单刚(局部坐标)右手

14、坐标系06*E*IL22*E*IL0-6*E*II24*E*IL;returnzbzh.m%形成第i单元的坐标转换矩阵函数T(6,6)%人口参数：单元号，单元信息，结点坐标%出口参数：坐标转换矩阵(整体向局部投影)functionT=zbzh(i,LNODS,COORD)NL=LNODS(i,l);NR=LNODS(i,2);%左结点号%右结点号dx=COORD(NR,l)-COORD(NL,l);%x坐标差dy=COORD(NR,2)-COORD(NL,2);%y坐标差1.=sqrt(dx2+dy2);%单元长度c=dxL;%cosa(与X轴夹角余弦)s=dyL;T=CSo-Sc0001000000000%sinareturn00;.00;.00;.s0;.c0;.0U；