基本平面刚架MATLAB程序.doc

% 平面刚架MATLAB程序% 2003.9.16 2007.2.28 2008.4.1 2009.10 2011.10 2013.9 2014.09%*% 变量说明% NPOIN NELEM NVFIX NFPOIN NFPRES% 总结点数,单元数, 约束个数, 受力结点数, 非结点力数% COORD LNODS YOUNG% 结构节点坐标数组, 单元定义数组, 弹性模量A% FPOIN FPRES FORCE FIXED% 结点力数组，非结点力数组，总体荷载向量, 约束信息数组% HK DISP% 总体刚度矩阵，结点位移向量%*format short e %设定输出类型clear %清除内存变量FP1=fopen(XX.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); % 张成总荷载向量并清零%形成总刚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); End % 单刚3 x 3子块叠加到总刚中 endend% 由结点力与非结点力生成总荷载向量列阵for i=1:NFPOIN % 对结点荷载个数进行循环 N1=FPOIN(i,1); % 作用荷载的结点号 N1=N1*3-3; % 该结点号对应第一个位移编码 - 1 for j=1:3 FORCE(N1+j)=FORCE(N1+j)+FPOIN(i,j+1);%取结点荷载 end end% 计算由非结点荷载引起的等效结点荷载for i=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% 总刚、总荷载进行边界条件处理for j=1:NVFIX % 对约束个数进行循环 N1=FIXED(j);HK(1:3*NPOIN,N1)=0; HK(N1,1:3*NPOIN)=0; HK(N1,N1)=1; % 将零位移约束对应的行、列变成零，主元变成1 FORCE(N1)=0;end%-DISP=HKFORCE % 方程求解，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:NFPRES