基本平面刚架MATLAB程序
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
% 平面刚架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子块叠加到总刚中
end
end
% 由结点力与非结点力生成总荷载向量列阵
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=HK\FORCE % 方程求解,HK先求逆再与力向量左乘
% 求结构各个单元内力
EDISP=zeros(6,1); % 单元位移列向量清零