2016基本平面刚架各种荷载MATLAB程序

% 平面刚架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子块叠加到总刚中

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,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;

% 将零位移约束对应的行、列变成零，主元变成1

FORCE(N1)=0;

end

%---------------------------------------------------------

DISP=HK\FORCE % 方程求解，HK先求逆再与力向量左乘