% 平面刚架MATLAB程序% 2003.9.16 2007.2.28 2008.4.1 2009.10 2011.10 2013.9 2014.09 2016.03%*************************************************% 变量说明% NPOIN NELEM NVFIX NFPOIN NFPRES% 总结点数,单元数, 约束个数, 受力结点数, 非结点力数% COORD LNODS YOUNG% 结构节点坐标数组, 单元定义数组, 弹性模量% FPOIN FPRES FORCE FIXED% 结点力数组，非结点力数组，总体荷载向量, 约束信息数组% HK DISP% 总体刚度矩阵，结点位移向量%**************************************************format short e %设定输出类型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',[6,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',[7,NFPRES])' % 均布力（共计% NFPRES 组）：单元号、荷载类型、荷载大小、距离左端长度，温差=（下端-上端）梯形上边。

下边(改)% 荷载类型1-均布荷载2-横向集中力3-纵向集中力4-三角形荷载5-温度荷载6-梯形荷载FIXED=fscan f(FP1,'%f',NVFIX)'% 约束信息：约束对应的位移编码（共计NVFIX 组）%---------------------------------------------------------HK=zeros(3*NPOIN,3*NPOIN); % 张成总刚矩阵并清零FORCE=zeros(3*NPOIN,1); % 张成总荷载向量并清零%形成总刚for i=1:NELEM % 对单元个数循环% 生成局部单刚(局部坐标) 右手坐标系EK=ele_EK(i,LNODS,COORD,YOUNG);T=zbzh(i,LNODS,COORD);% 坐标转换矩阵EKT=T’*EK*T；% 生成整体单刚(整体坐标系)% 组成总刚按3*3子块加入总刚中（共计4块)for j=1:2 %对行进行循环---按结点号循环N1=LNODS(i,j)*3; % j结点第3个位移的整体编码for k=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);% 单刚3 x 3子块叠加到总刚中endend% 由结点力与非结点力生成总荷载向量列阵for i=1:NFPOIN % 对结点荷载个数进行循环N1=FPOIN(i,1); % 作用荷载的结点号N1=N1*3-3; % 该结点号对应第一个位移编码- 1for j=1:3FORCE(N1+j)=FORCE(N1+j)+FPOIN(i,j+1);%取结点荷载endend% 计算由非结点荷载引起的等效结点荷载for i=1:NFPRES % 对非结点荷载个数进行循环F0=ele_FPRES(i,FPRES,LNODS,COORD,YOUNG); %计算单元固端力% 对单元局部杆端力要进行坐标转换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% 总刚、总荷载进行边界条件处理for j=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); % 单元位移列向量清零for i=1:NELEM % 对单元个数进行循环for j=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; %计算局部坐标杆端力（由结点位移产生）for j=1:NFPRESif FPRES(j,1) == i %成立时，当前单元上有非结点荷载F0=ele_FPRES(j,FPRES,LNODS,COORD,YOUNG);%单元固端力FE=FE+F0; % 考虑由非结点荷载引起的杆端力endendFE % 打印杆端力endend%-------------------------------------------------------ele_FPRES.m %计算单元固端力函数（正方向：X向右Y向上M逆时针）% 入口参数：荷载序号，荷载信息，单元信息，结点坐标% 出口参数：单元固端力——左右两端的轴力、剪力、弯矩function F0=ele_FPRES(iFPRES,FPRES,LNODS,COORD, E)ele=FPRES(iFPRES,1); %取荷载所在的单元号G=FPRES(iFPRES,3); %单元荷载大小C=FPRES(iFPRES,4); %单元荷载与左端距离W= FPRES(iFPRES,5); %单元下上端温差S= FPRES(iFPRES,6); %梯形长边荷载X= FPRES(iFPRES,6); %梯形短边荷载H= LNODS(ele,6);%单元截面高度P = LNODS(ele,5);%单元截面线膨胀系数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]; %单元固端力清零switch FPRES(iFPRES,2)case 1 %均布荷载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;case 2 %横向集中力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;case 3 %纵向集中力F0(1)=-G*B1;F0(4)=-G*C1;case 4 %三角形荷载F0(2)=-7*G*L/20;F0(3)=G*L^2/20;F0(5)=3*G*L/20;F0(6)=G*L*L/30;case 5 %温度荷载F0(3)=E*I*W*P/H;F0(6)=-E*I*W*P/H;case 6 %梯形荷载（改）F1(2)=-X*C*(2-2*C2+C3)/2.0;F1(3)=-X*C*C*(6-8*C1+3*C2)/12.0;F2(5)=-X*C-F0(2);F2(6)=X*C*C*C1*(4-3*C1)/12.0;F3(2)=-7*(S-G)*L/20;F3(3)=(S-G)*L^2/20;F4(5)=3*(S-G)*L/20;F4(6)=(S-G)*L*L/30;F0(2)= F1(2)+ F3(2);F0(3)= F1(3)+ F3(3);F0(5)= F2(5)+ F4(5);F0(6)= F2(6)+ F4(6);endreturnele_EK.m % 计算单元刚度矩阵函数EK% 入口参数：单元号、单元信息数组、结点坐标、弹性模量% 出口参数：局部单元刚度矩阵EKfunction EK=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/L 0 0 -E*A/L 0 0;...0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 6*E*I/L^2;...0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L;...-E*A/L 0 0 E*A/L 0 0;...0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -6*E*I/L^2;...0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 4*E*I/L]; return%--------------------------------------------------------zbzh.m % 形成第i单元的坐标转换矩阵函数T(6,6)% 入口参数：单元号，单元信息，结点坐标% 出口参数：坐标转换矩阵（整体向局部投影）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); % 单元长度c=dx/L; % cos a (与x 轴夹角余弦）s=dy/L; % sin aT=[ c s 0 0 0 0;...-s c 0 0 0 0;...0 0 1 0 0 0;...0 0 0 c s 0;...0 0 0 -s c 0;...0 0 0 0 0 1];return。